CR2RE Pipeline Reference Manual 1.7.0
hdrl_resample-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#include <cpl.h>
29#include <stdint.h>
30#include <math.h>
31#include <assert.h>
32#include <hdrl_utils.h>
33#include <hdrl_resample.h>
34#include <hdrl_resample.c>
35#include <hdrl_random.h>
36#ifndef _OPENMP
37#define omp_get_max_threads() 1
38#else
39#include <omp.h>
40#endif
41
42/* The following code contains several unit tests to cover the hdrl_resample.c
43 * library. As we would like that any function has its own unit test for easy
44 * we used a naming convention as follow:
45 *
46 * 'test_function()' tests the function 'function()';
47 *
48 * There are a few exceptions: "in a few cases it was more convenient to have
49 * a single function testing several library functions
50 * Moreover we have some additional unit tests, testing special functionality
51 * which follows a different naming convention.
52 *
53 */
54/* -------------------------------- DEFINES --------------------------------- */
55#define RECIPE_NAME "hdrldemo_resample"
56
57/*Fine-grained accuracies*/
58#define HDRL_DELTA_COMPARE_VALUE_ABS CPL_MAX(HDRL_EPS_DATA, HDRL_EPS_ERROR) * 4.
59
60#define HDRL_EPS_TEST HDRL_EPS_DATA
61/* For resampling method definition */
62#define DRIZZLE_DOWN_SCALING_FACTOR_X 0.8
63#define DRIZZLE_DOWN_SCALING_FACTOR_Y 0.8
64#define DRIZZLE_DOWN_SCALING_FACTOR_Z 0.8
65#define RENKA_CRITICAL_RADIUS 1.25
66#define LANCZOS_KERNEL_SIZE 2
67#define LOOP_DISTANCE 1.0
68/* Image value */
69#define HDRL_FLUX_ADU 100
70/* For WCS definition From SINFONI cube example 1:
71 * SINFO.2005-08-22T07:47:54.305.fits */
72#define HDRL_SCALE_Z 500
73#define HDRL_CD11 -3.47222e-05
74#define HDRL_CD12 0.0
75#define HDRL_CD21 0.0
76#define HDRL_CD22 3.47222e-05
77#define HDRL_CD13 0.0
78#define HDRL_CD31 0.0
79#define HDRL_CD23 0.0
80#define HDRL_CD32 0.0
81#define HDRL_CD33 2.45e-10*HDRL_SCALE_Z
82
83#define HDRL_CDELT1 fabs(HDRL_CD11)
84#define HDRL_CDELT2 fabs(HDRL_CD22)
85#define HDRL_CDELT3 fabs(HDRL_CD33)
86
87#define HDRL_CRPIX1 33.5
88#define HDRL_CRPIX2 33.5
89#define HDRL_CRPIX3 1.
90
91#define HDRL_CRVAL1 48.0706
92#define HDRL_CRVAL2 -20.6219
93#define HDRL_CRVAL3 1.9283e-06
94
95#define HDRL_RA 48.070
96#define HDRL_DEC -20.621
97#define HDRL_RA_MIN 48.069416667
98#define HDRL_RA_MAX 48.0718125
99#define HDRL_DEC_MIN -20.6229925
100#define HDRL_DEC_MAX -20.620708611
101
102#define HDRL_LAMBDA_MIN 1.9283e-06
103#define HDRL_LAMBDA_MAX 2.47146e-06
104
105/* Image sizes */
106#define HDRL_SIZE_X 50
107#define HDRL_SIZE_Y 50
108#define HDRL_SIZE_Z 3
109
110
111/*----------------------------------------------------------------------------*/
116/*----------------------------------------------------------------------------*/
117static hdrl_image*
118hdrl_resample_util_hdrl_image_create(cpl_image* data, cpl_image* error,
119 cpl_image* bpm){
120 if(bpm != NULL) {
121 cpl_mask* mask = cpl_mask_threshold_image_create(bpm, 0, INT_MAX);
122 cpl_image_reject_from_mask(data, mask);
123 cpl_mask_delete(mask);
124 }
125 hdrl_image* hima = hdrl_image_create(data, error);
126
127 return hima;
128}
129
130/*----------------------------------------------------------------------------*/
135/*----------------------------------------------------------------------------*/
136static cpl_error_code test_hdrl_resample_compute2D_multiple(void)
137{
138
139 /* Data background always 1 and inner 9 pixel 49 - for error tracing*/
140
141 cpl_image * img_data = cpl_image_new(9, 9, CPL_TYPE_DOUBLE);
142 cpl_image * img_bpm = cpl_image_new(9, 9, CPL_TYPE_INT);
143
144 /* Set all data invalid */
145 cpl_image_add_scalar(img_bpm, 1.);
146
147 /* Set data in the center of the image */
148 cpl_image_set(img_data, 4, 4, 48.);
149 cpl_image_set(img_data, 5, 4, 48.);
150 cpl_image_set(img_data, 6, 4, 48.);
151 cpl_image_set(img_data, 4, 5, 48.);
152 cpl_image_set(img_data, 5, 5, 48.);
153 cpl_image_set(img_data, 6, 5, 48.);
154 cpl_image_set(img_data, 4, 6, 48.);
155 cpl_image_set(img_data, 5, 6, 48.);
156 cpl_image_set(img_data, 6, 6, 48.);
157 /* Adding 1 and creating the errors */
158 cpl_image_add_scalar(img_data, 1.0);
159 cpl_image * img_error=cpl_image_power_create(img_data, 0.5);
160
161 /* Set data in the center as valid */
162 cpl_image_set(img_bpm, 4, 4, 0);
163 cpl_image_set(img_bpm, 5, 4, 0);
164 cpl_image_set(img_bpm, 6, 4, 0);
165 cpl_image_set(img_bpm, 4, 5, 0);
166 cpl_image_set(img_bpm, 5, 5, 0);
167 cpl_image_set(img_bpm, 6, 5, 0);
168 cpl_image_set(img_bpm, 4, 6, 0);
169 cpl_image_set(img_bpm, 5, 6, 0);
170 cpl_image_set(img_bpm, 6, 6, 0);
171
172 /*Build the header for the wcs */
173 cpl_propertylist* plist = cpl_propertylist_new();
174
175 cpl_propertylist_append_int(plist, "NAXIS" , 2);
176 cpl_propertylist_append_int(plist, "NAXIS1" , 9);
177 cpl_propertylist_append_int(plist, "NAXIS2" , 9);
178 cpl_propertylist_append_double(plist, "CD1_1" , -0.01 );
179 cpl_propertylist_append_double(plist, "CD1_2" , 0. );
180 cpl_propertylist_append_double(plist, "CD2_1" , 0. );
181 cpl_propertylist_append_double(plist, "CD2_2" , 0.01 );
182 cpl_propertylist_append_double(plist, "CRPIX1" , 4.5 );
183 cpl_propertylist_append_double(plist, "CRPIX2" , 4.5 );
184 cpl_propertylist_append_double(plist, "CRVAL1" , 359.8 );
185 cpl_propertylist_append_double(plist, "CRVAL2" , 10.0 );
186 cpl_propertylist_append_string(plist, "CTYPE1" , "RA---TAN");
187 cpl_propertylist_append_string(plist, "CTYPE2" , "DEC--TAN");
188 cpl_propertylist_append_string(plist, "CUNIT1" , "deg" );
189 cpl_propertylist_append_string(plist, "CUNIT2" , "deg" );
190
191 hdrl_random_state * rastate = hdrl_random_state_new(1, NULL);
192 hdrl_random_state * decstate = hdrl_random_state_new(1, NULL);
193
194 /* Please be aware on the different accuracies set for the errors for the
195 * different methods */
196
197 /* generate hdrl image */
198 hdrl_image*
199 hima = hdrl_resample_util_hdrl_image_create(img_data, img_error, img_bpm);
200 const cpl_size iterations = 500;
201 for (cpl_size i = 0; i < iterations; i++) {
202 double rarandom = hdrl_random_uniform_double(rastate, 0., 360.);
203 double decrandom = hdrl_random_uniform_double(decstate, -89., 89.);
204 cpl_msg_info(cpl_func,"ra-random: %g, dec-random: %g ", rarandom,
205 decrandom);
206 cpl_propertylist_update_double(plist, "CRVAL1" , rarandom);
207 cpl_propertylist_update_double(plist, "CRVAL2" , decrandom);
208 cpl_wcs * wcs = cpl_wcs_new_from_propertylist(plist);
209
210 /* Construct the table */
211 cpl_table * table = hdrl_resample_image_to_table(hima, wcs);
212
213 hdrl_parameter *aParams_outputgrid = NULL;
214
215 /* Define the output grid */
216 aParams_outputgrid = hdrl_resample_parameter_create_outgrid2D(0.01, 0.01);
217 hdrl_resample_result *result = NULL;
218 hdrl_parameter *aParams_method = NULL;
219 int rej;
220
221 /* ---------------------------lanczos------------------------------- */
222
223 /* Define the method */
224 aParams_method = hdrl_resample_parameter_create_lanczos(1, CPL_FALSE, 2);
225 /* Do the resampling */
226 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
227 cpl_test_error(CPL_ERROR_NONE);
228 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 0)),
229 5, 5, &rej), 49., 1e-6);
230 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 0)),
231 5, 5, &rej), 7., 0.05);
232 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
233 5, 5), 0);
234 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
235 1, 1), 1);
236 hdrl_parameter_delete(aParams_method);
237 hdrl_resample_result_delete(result);
238
239 /* -----------------------------drizzle----------------------------- */
240
241 /* Define the method */
242 aParams_method = hdrl_resample_parameter_create_drizzle(1, CPL_FALSE,
243 0.8, 0.8, 0.8);
244 /* Do the resampling */
245
246 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
247 cpl_test_error(CPL_ERROR_NONE);
248 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 0)),
249 5, 5, &rej), 49., 1e-6);
250 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 0)),
251 5, 5, &rej), 7., 0.001);
252 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
253 5, 5), 0);
254 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
255 1, 1), 1);
256
257 hdrl_parameter_delete(aParams_method);
258 hdrl_resample_result_delete(result);
259
260 /* -----------------------------linear------------------------------ */
261
262 /* Define the method */
263 aParams_method = hdrl_resample_parameter_create_linear(1, CPL_FALSE);
264 /* Do the resampling */
265
266 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
267 cpl_test_error(CPL_ERROR_NONE);
268 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 0)),
269 5, 5, &rej), 49., 1e-6);
270 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 0)),
271 5, 5, &rej), 7., 0.3);
272 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
273 5, 5), 0);
274 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
275 1, 1), 1);
276
277 hdrl_parameter_delete(aParams_method);
278 hdrl_resample_result_delete(result);
279
280 /* -----------------------------quadratic--------------------------- */
281
282 /* Define the method */
283 aParams_method = hdrl_resample_parameter_create_quadratic(1, CPL_FALSE);
284 /* Do the resampling */
285
286 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
287 cpl_test_error(CPL_ERROR_NONE);
288 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 0)),
289 5, 5, &rej), 49., 1e-6);
290 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 0)),
291 5, 5, &rej), 7., 0.02);
292 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
293 5, 5), 0);
294 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
295 1, 1), 1);
296
297 hdrl_parameter_delete(aParams_method);
298 hdrl_resample_result_delete(result);
299
300 /* -----------------------------renka------------------------------- */
301
302 /* Define the method */
303 aParams_method = hdrl_resample_parameter_create_renka(1, CPL_FALSE,
304 1.25);
305 /* Do the resampling */
306
307 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
308 cpl_test_error(CPL_ERROR_NONE);
309 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 0)),
310 5, 5, &rej), 49., 1e-6);
311 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 0)),
312 5, 5, &rej), 7., 0.01);
313 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
314 5, 5), 0);
315 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
316 1, 1), 1);
317
318 hdrl_parameter_delete(aParams_method);
319 hdrl_resample_result_delete(result);
320
321 /* -----------------------------nearest----------------------------- */
322
323 /* Define the method */
324 aParams_method = hdrl_resample_parameter_create_nearest();
325 /* Do the resampling */
326 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
327 cpl_test_error(CPL_ERROR_NONE);
328 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 0)),
329 5, 5, &rej), 49., 1e-6);
330 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 0)),
331 5, 5, &rej), 7., 0.00001);
332 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
333 5, 5), 0);
334 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
335 1, 1), 1);
336
337 hdrl_parameter_delete(aParams_method);
338 hdrl_resample_result_delete(result);
339
340 /* ----------------------------------------------------------------- */
341
342 hdrl_parameter_delete(aParams_outputgrid);
343 cpl_table_delete(table);
344 cpl_wcs_delete(wcs);
345 } /* Loop ends here */
346 hdrl_image_delete(hima);
347
348 /* ---------------- For manual updates and tests ----------------------- */
349
350 /* Test image spanning over ra 360 degree */
351 cpl_propertylist_update_double(plist, "CRVAL1" , 0.01);
352 cpl_propertylist_update_double(plist, "CRVAL2" , 20.1);
353
354 /* Test image rotated by 45 degree */
355 cpl_propertylist_update_double(plist, "CD1_1" , cos(45) );
356 cpl_propertylist_update_double(plist, "CD1_2" , -sin(45) );
357 cpl_propertylist_update_double(plist, "CD2_1" , sin(45) );
358 cpl_propertylist_update_double(plist, "CD2_2" , cos(45) );
359
360
361 cpl_wcs * wcs = cpl_wcs_new_from_propertylist(plist);
362
363 /* Construct the table */
364 hima = hdrl_resample_util_hdrl_image_create(img_data, img_error, img_bpm);
365
366 cpl_table * table = hdrl_resample_image_to_table(hima, wcs);
367
368 hdrl_image_delete(hima);
369 hdrl_parameter *aParams_outputgrid = NULL;
370
371 /* Define the output grid */
372 aParams_outputgrid = hdrl_resample_parameter_create_outgrid2D(cos(45) * 3., cos(45) * 3.);
373
374 hdrl_resample_result *result = NULL;
375
376 hdrl_parameter *aParams_method = NULL;
377 int rej;
378
379 /* Define the method */
380 aParams_method = hdrl_resample_parameter_create_lanczos(2, CPL_FALSE, 2);
381 /* Do the resampling */
382 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
383 cpl_test_error(CPL_ERROR_NONE);
384 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 0)),
385 5, 5, &rej), 49., 1e-6);
386 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 0)),
387 5, 5, &rej), 7., 0.3);
388 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
389 5, 5), 0);
390 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 0)),
391 1, 1), 1);
392
393 /*Save the original and resampled image/error/bpm*/
394 cpl_image_save(img_data, "image_orig.fits", CPL_TYPE_FLOAT, plist,
395 CPL_IO_CREATE);
396
397 cpl_image_save(img_error, "image_orig.fits", CPL_TYPE_FLOAT, plist,
398 CPL_IO_EXTEND);
399
400 cpl_image_save(img_bpm, "image_orig.fits", CPL_TYPE_INT, plist,
401 CPL_IO_EXTEND);
402
403 cpl_image_save(hdrl_image_get_image(hdrl_imagelist_get(result->himlist, 0)), "image_resampled.fits",
404 CPL_TYPE_FLOAT, result->header, CPL_IO_CREATE);
405 cpl_image_save(hdrl_image_get_error(hdrl_imagelist_get(result->himlist, 0)), "image_resampled.fits",
406 CPL_TYPE_FLOAT, result->header, CPL_IO_EXTEND);
407 cpl_mask_save(hdrl_image_get_mask(hdrl_imagelist_get(result->himlist, 0)), "image_resampled.fits",
408 result->header, CPL_IO_EXTEND);
409
410 /*Print the wcs from the original and resampled cube*/
411 hdrl_resample_wcs_print(wcs);
412
413 cpl_wcs_delete(wcs);
414 wcs = cpl_wcs_new_from_propertylist(result->header);
415 hdrl_resample_wcs_print(wcs);
416
417 hdrl_parameter_delete(aParams_method);
418 hdrl_resample_result_delete(result);
419 hdrl_parameter_delete(aParams_outputgrid);
420 cpl_table_delete(table);
421 cpl_wcs_delete(wcs);
422
423 /* --------------------------------------------------------------------- */
424
425 /* free the remaining memory */
426 cpl_propertylist_delete(plist);
427 cpl_image_delete(img_data);
428 cpl_image_delete(img_error);
429 cpl_image_delete(img_bpm);
430 hdrl_random_state_delete(rastate);
431 hdrl_random_state_delete(decstate);
432
433 cpl_test_error(CPL_ERROR_NONE);
434 return cpl_error_get_code();
435}
436
437static hdrl_imagelist*
438hdrl_resample_util_hdrl_imagelist_create(cpl_imagelist* dlist,
439 cpl_imagelist* elist, cpl_imagelist* qlist) {
440
441 cpl_ensure(dlist != NULL, CPL_ERROR_NULL_INPUT, NULL);
442 cpl_size size = cpl_imagelist_get_size(dlist);
443 if(qlist != NULL) {
444 for(cpl_size k = 0; k < size; k++) {
445 cpl_image* data = cpl_imagelist_get(dlist, k);
446
447 cpl_image* qual = cpl_imagelist_get(qlist, k);
448
449 /*we use INT_MAX instead of 1.1 as some pipeline
450 * may use pixel codes as qualifier */
451 cpl_mask* mask = cpl_mask_threshold_image_create(qual, 0, INT_MAX);
452
453 cpl_image_reject_from_mask(data, mask);
454 cpl_mask_delete(mask);
455 cpl_imagelist_set(dlist, data, k);
456 }
457 }
458 hdrl_imagelist* hlist = hdrl_imagelist_create(dlist, elist);
459
460 return hlist;
461
462}
463/*----------------------------------------------------------------------------*/
468/*----------------------------------------------------------------------------*/
469static cpl_error_code test_hdrl_resample_compute3D_multiple(void)
470{
471 cpl_test_error(CPL_ERROR_NONE);
472 cpl_imagelist * imglist_data = cpl_imagelist_new();
473 cpl_imagelist * imglist_error = cpl_imagelist_new();
474 cpl_imagelist * imglist_bpm = cpl_imagelist_new();
475 cpl_test_error(CPL_ERROR_NONE);
476
477 /* Data background always 0 and inner 9 pixel ramp up from 1 to 50 and then
478 * back to 1 */
479 cpl_test_error(CPL_ERROR_NONE);
480 for (cpl_size var = 0; var <99 ; ++var) {
481 cpl_image * data = cpl_image_new(9, 9, CPL_TYPE_DOUBLE);
482 cpl_image * bpm = cpl_image_new(9, 9, CPL_TYPE_INT);
483 /* Set all data invalid */
484 cpl_image_add_scalar(bpm, 1.);
485
486 if(var < 50) {
487 cpl_image_set(data, 4, 4, var + 1.);
488 cpl_image_set(data, 5, 4, var + 1.);
489 cpl_image_set(data, 6, 4, var + 1.);
490 cpl_image_set(data, 4, 5, var + 1.);
491 cpl_image_set(data, 5, 5, var + 1.);
492 cpl_image_set(data, 6, 5, var + 1.);
493 cpl_image_set(data, 4, 6, var + 1.);
494 cpl_image_set(data, 5, 6, var + 1.);
495 cpl_image_set(data, 6, 6, var + 1.);
496 cpl_image * errors=cpl_image_power_create(data, 0.5);
497
498 /* Set data in the center as valid */
499 cpl_image_set(bpm, 4, 4, 0);
500 cpl_image_set(bpm, 5, 4, 0);
501 cpl_image_set(bpm, 6, 4, 0);
502 cpl_image_set(bpm, 4, 5, 0);
503 cpl_image_set(bpm, 5, 5, 0);
504 cpl_image_set(bpm, 6, 5, 0);
505 cpl_image_set(bpm, 4, 6, 0);
506 cpl_image_set(bpm, 5, 6, 0);
507 cpl_image_set(bpm, 6, 6, 0);
508
509 cpl_imagelist_set(imglist_data, data, var);
510 cpl_imagelist_set(imglist_error, errors, var);
511 cpl_imagelist_set(imglist_bpm, bpm, var);
512 }
513 if(var >= 50) {
514 cpl_image_set(data, 4, 4, 99. - var);
515 cpl_image_set(data, 5, 4, 99. - var);
516 cpl_image_set(data, 6, 4, 99. - var);
517 cpl_image_set(data, 4, 5, 99. - var);
518 cpl_image_set(data, 5, 5, 99. - var);
519 cpl_image_set(data, 6, 5, 99. - var);
520 cpl_image_set(data, 4, 6, 99. - var);
521 cpl_image_set(data, 5, 6, 99. - var);
522 cpl_image_set(data, 6, 6, 99. - var);
523 cpl_image * errors=cpl_image_power_create(data, 0.5);
524
525 /* Set data in the center as valid */
526 cpl_image_set(bpm, 4, 4, 0);
527 cpl_image_set(bpm, 5, 4, 0);
528 cpl_image_set(bpm, 6, 4, 0);
529 cpl_image_set(bpm, 4, 5, 0);
530 cpl_image_set(bpm, 5, 5, 0);
531 cpl_image_set(bpm, 6, 5, 0);
532 cpl_image_set(bpm, 4, 6, 0);
533 cpl_image_set(bpm, 5, 6, 0);
534 cpl_image_set(bpm, 6, 6, 0);
535
536 cpl_imagelist_set(imglist_data, data, var);
537 cpl_imagelist_set(imglist_error, errors, var);
538 cpl_imagelist_set(imglist_bpm, bpm, var);
539 }
540 }
541 cpl_test_error(CPL_ERROR_NONE);
542 /*Build the header for the wcs */
543
544 cpl_propertylist* plist = cpl_propertylist_new();
545
546 cpl_propertylist_append_int(plist, "NAXIS" , 3 );
547 cpl_propertylist_append_int(plist, "NAXIS1" , 9 );
548 cpl_propertylist_append_int(plist, "NAXIS2" , 9 );
549 cpl_propertylist_append_int(plist, "NAXIS3" , 99 );
550 cpl_propertylist_append_double(plist, "CD1_1" , -0.01 );
551 cpl_propertylist_append_double(plist, "CD1_2" , 0. );
552 cpl_propertylist_append_double(plist, "CD2_1" , 0. );
553 cpl_propertylist_append_double(plist, "CD2_2" , 0.01 );
554 cpl_propertylist_append_double(plist, "CRPIX1" , 4.5 );
555 cpl_propertylist_append_double(plist, "CRPIX2" , 4.5 );
556 cpl_propertylist_append_double(plist, "CRVAL1" , 48.0 );
557 cpl_propertylist_append_double(plist, "CRVAL2" , -20.0 );
558 cpl_propertylist_append_string(plist, "CTYPE1" , "RA---TAN");
559 cpl_propertylist_append_string(plist, "CTYPE2" , "DEC--TAN");
560 cpl_propertylist_append_string(plist, "CUNIT1" , "deg" );
561 cpl_propertylist_append_string(plist, "CUNIT2" , "deg" );
562 cpl_propertylist_append_double(plist, "CD1_3" , 0. );
563 cpl_propertylist_append_double(plist, "CD2_3" , 0. );
564 cpl_propertylist_append_double(plist, "CD3_1" , 0. );
565 cpl_propertylist_append_double(plist, "CD3_2" , 0. );
566 cpl_propertylist_append_double(plist, "CD3_3" , 1.0 );
567 cpl_propertylist_append_double(plist, "CRPIX3" , 1.0 );
568 cpl_propertylist_append_double(plist, "CRVAL3" , 1.0 );
569 cpl_propertylist_append_string(plist, "CTYPE3" , "WAVE" );
570 cpl_propertylist_append_string(plist, "CUNIT3" , "m" );
571 cpl_test_error(CPL_ERROR_NONE);
572
573 hdrl_random_state * rastate = hdrl_random_state_new(1, NULL);
574 hdrl_random_state * decstate = hdrl_random_state_new(1, NULL);
575
576 /* Please be aware on the different accuracies set for the errors for the
577 * different methods */
578 hdrl_imagelist* hlist =
579 hdrl_resample_util_hdrl_imagelist_create(imglist_data,
580 imglist_error, imglist_bpm);
581 const cpl_size iterations = 20;
582 for (cpl_size i = 0; i < iterations; i++) {
583 double rarandom = hdrl_random_uniform_double(rastate, 0., 360.);
584 double decrandom = hdrl_random_uniform_double(decstate, -89., 89.);
585 cpl_msg_info(cpl_func,"ra-random: %g, dec-random: %g ", rarandom,
586 decrandom);
587 cpl_propertylist_update_double(plist, "CRVAL1" , rarandom);
588 cpl_propertylist_update_double(plist, "CRVAL2" , decrandom);
589 cpl_wcs * wcs = cpl_wcs_new_from_propertylist(plist);
590 cpl_test_error(CPL_ERROR_NONE);
591 /* Construct the table */
592 cpl_table * table = hdrl_resample_imagelist_to_table(hlist, wcs);
593
594 cpl_test_error(CPL_ERROR_NONE);
595 hdrl_parameter *aParams_outputgrid = NULL;
596 /* Define the output grid */
597 aParams_outputgrid =
598 hdrl_resample_parameter_create_outgrid3D(0.01, 0.01, 1);
599 cpl_test_nonnull(aParams_outputgrid);
600 cpl_test_error(CPL_ERROR_NONE);
601
602 hdrl_resample_result *result = NULL;
603 hdrl_parameter *aParams_method = NULL;
604 int rej;
605
606 /* ---------------------------lanczos------------------------------- */
607
608 /* Define the method */
609 aParams_method = hdrl_resample_parameter_create_lanczos(1, CPL_FALSE, 2);
610 /* Do the resampling */
611 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
612 cpl_test_error(CPL_ERROR_NONE);
613 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 48)),
614 5, 5, &rej), 49., 1e-6);
615 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 48)),
616 5, 5, &rej), 7., 0.05);
617 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
618 5, 5), 0);
619 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
620 1, 1), 1);
621 hdrl_parameter_delete(aParams_method);
622 hdrl_resample_result_delete(result);
623
624
625 /* -----------------------------drizzle----------------------------- */
626
627 /* Define the method */
628 aParams_method = hdrl_resample_parameter_create_drizzle(1, CPL_FALSE,
629 0.8, 0.8, 0.8);
630 /* Do the resampling */
631
632 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
633 cpl_test_error(CPL_ERROR_NONE);
634 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 48)),
635 5, 5, &rej), 49., 1e-6);
636 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 48)),
637 5, 5, &rej), 7., 0.001);
638 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
639 5, 5), 0);
640 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
641 1, 1), 1);
642
643 hdrl_parameter_delete(aParams_method);
644 hdrl_resample_result_delete(result);
645
646 /* -----------------------------linear------------------------------ */
647
648 /* Define the method */
649 aParams_method = hdrl_resample_parameter_create_linear(1, CPL_FALSE);
650 /* Do the resampling */
651
652 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
653 cpl_test_error(CPL_ERROR_NONE);
654 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 48)),
655 5, 5, &rej), 49., 1e-6);
656 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 48)),
657 5, 5, &rej), 7., 0.3);
658 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
659 5, 5), 0);
660 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
661 1, 1), 1);
662
663 hdrl_parameter_delete(aParams_method);
664 hdrl_resample_result_delete(result);
665
666 /* -----------------------------quadratic--------------------------- */
667
668 /* Define the method */
669 aParams_method = hdrl_resample_parameter_create_quadratic(1, CPL_FALSE);
670 /* Do the resampling */
671
672 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
673 cpl_test_error(CPL_ERROR_NONE);
674 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 48)),
675 5, 5, &rej), 49., 1e-6);
676 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 48)),
677 5, 5, &rej), 7., 0.02);
678 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
679 5, 5), 0);
680 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
681 1, 1), 1);
682
683 hdrl_parameter_delete(aParams_method);
684 hdrl_resample_result_delete(result);
685
686 /* -----------------------------renka------------------------------- */
687
688 /* Define the method */
689 aParams_method = hdrl_resample_parameter_create_renka(1, CPL_FALSE,
690 1.25);
691 /* Do the resampling */
692
693 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
694 cpl_test_error(CPL_ERROR_NONE);
695 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 48)),
696 5, 5, &rej), 49., 1e-6);
697 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 48)),
698 5, 5, &rej), 7., 0.01);
699 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
700 5, 5), 0);
701 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
702 1, 1), 1);
703
704 hdrl_parameter_delete(aParams_method);
705 hdrl_resample_result_delete(result);
706
707 /* -----------------------------nearest----------------------------- */
708
709 /* Define the method */
710 aParams_method = hdrl_resample_parameter_create_nearest();
711 /* Do the resampling */
712 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
713 cpl_test_error(CPL_ERROR_NONE);
714 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 48)),
715 5, 5, &rej), 49., 1e-6);
716 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 48)),
717 5, 5, &rej), 7., 0.00001);
718 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
719 5, 5), 0);
720 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
721 1, 1), 1);
722
723 hdrl_parameter_delete(aParams_method);
724 hdrl_resample_result_delete(result);
725
726 /* ----------------------------------------------------------------- */
727
728 hdrl_parameter_delete(aParams_outputgrid);
729 cpl_table_delete(table);
730 cpl_wcs_delete(wcs);
731 } /* Loop ends here */
732
733
734 /* ---------------- For manual updates and tests ----------------------- */
735
736 /* Test image spanning over ra 360 degree */
737 cpl_propertylist_update_double(plist, "CRVAL1" , 0.03);
738 cpl_propertylist_update_double(plist, "CRVAL2" , 0.1);
739 cpl_wcs * wcs = cpl_wcs_new_from_propertylist(plist);
740
741 /* Construct the table */
742 cpl_table * table = hdrl_resample_imagelist_to_table(hlist, wcs);
743 hdrl_imagelist_delete(hlist);
744 hdrl_parameter *aParams_outputgrid = NULL;
745 /* Define the output grid */
746 aParams_outputgrid =
747 hdrl_resample_parameter_create_outgrid3D(0.01, 0.01, 1);
748 cpl_test_nonnull(aParams_outputgrid);
749 cpl_test_error(CPL_ERROR_NONE);
750
751 hdrl_resample_result *result = NULL;
752 hdrl_parameter *aParams_method = NULL;
753 int rej;
754
755 /* Define the method */
756 aParams_method = hdrl_resample_parameter_create_lanczos(1, CPL_FALSE, 2);
757 /* Do the resampling */
758 result = hdrl_resample_compute(table, aParams_method, aParams_outputgrid, wcs);
759 cpl_test_error(CPL_ERROR_NONE);
760 cpl_test_rel(cpl_image_get(hdrl_image_get_image_const(hdrl_imagelist_get_const(result->himlist, 48)),
761 5, 5, &rej), 49., 1e-6);
762 cpl_test_rel(cpl_image_get(hdrl_image_get_error_const(hdrl_imagelist_get_const(result->himlist, 48)),
763 5, 5, &rej), 7., 0.05);
764 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
765 5, 5), 0);
766 cpl_test_eq(cpl_mask_get(hdrl_image_get_mask_const(hdrl_imagelist_get_const(result->himlist, 48)),
767 1, 1), 1);
768
769 cpl_test_error(CPL_ERROR_NONE);
770 /*Save the original cube*/
771 cpl_propertylist_save(plist, "cube_orig.fits", CPL_IO_CREATE);
772 cpl_imagelist_save(imglist_data, "cube_orig.fits", CPL_TYPE_FLOAT, plist,
773 CPL_IO_EXTEND);
774 cpl_imagelist_save(imglist_bpm, "cube_orig.fits", CPL_TYPE_INT, plist,
775 CPL_IO_EXTEND);
776 cpl_imagelist_save(imglist_error, "cube_orig.fits", CPL_TYPE_FLOAT,
777 plist, CPL_IO_EXTEND);
778
779 /*Save the resampled cube */
780
781 cpl_propertylist_save(plist, "cube_resampled.fits", CPL_IO_CREATE);
782 cpl_test_error(CPL_ERROR_NONE);
783
784 cpl_imagelist * ilistdata = cpl_imagelist_new();
785 cpl_imagelist * ilisterrors = cpl_imagelist_new();
786 cpl_imagelist * ilistbpm = cpl_imagelist_new();
787
788 cpl_size planes = hdrl_imagelist_get_size(result->himlist);
789 for (cpl_size i = 0; i < planes; i++) {
790 /* The images are duplicated to avoid errors reported by cppcheck */
791 cpl_imagelist_set(ilistdata, cpl_image_duplicate(hdrl_image_get_image(hdrl_imagelist_get(result->himlist, i))), i);
792 cpl_imagelist_set(ilisterrors, cpl_image_duplicate(hdrl_image_get_error(hdrl_imagelist_get(result->himlist, i))), i);
793 cpl_image * bpm = cpl_image_new_from_mask(hdrl_image_get_mask(hdrl_imagelist_get(result->himlist, i)));
794 cpl_imagelist_set(ilistbpm, bpm, i);
795 }
796
797 /* Save the data/bpm/errors in the right extension as a cube */
798 cpl_imagelist_save(ilistdata, "cube_resampled.fits", CPL_TYPE_FLOAT,
799 result->header, CPL_IO_EXTEND);
800 cpl_imagelist_save(ilistbpm, "cube_resampled.fits", CPL_TYPE_INT,
801 result->header, CPL_IO_EXTEND);
802 cpl_imagelist_save(ilisterrors, "cube_resampled.fits", CPL_TYPE_FLOAT,
803 result->header, CPL_IO_EXTEND);
804
805
806 cpl_imagelist_delete(ilistdata);
807 cpl_imagelist_delete(ilisterrors);
808 cpl_imagelist_delete(ilistbpm);
809
810 /*Print the wcs from the original and resampled cube*/
811 hdrl_resample_wcs_print(wcs);
812 cpl_wcs_delete(wcs);
813 wcs = cpl_wcs_new_from_propertylist(result->header);
814 hdrl_resample_wcs_print(wcs);
815
816 /* free the memory */
817 cpl_propertylist_delete(plist);
818 cpl_imagelist_delete(imglist_data);
819 cpl_imagelist_delete(imglist_error);
820 cpl_imagelist_delete(imglist_bpm);
821 hdrl_random_state_delete(rastate);
822 hdrl_random_state_delete(decstate);
823
824 hdrl_parameter_delete(aParams_method);
825 hdrl_parameter_delete(aParams_outputgrid);
826 hdrl_resample_result_delete(result);
827 cpl_table_delete(table);
828 cpl_wcs_delete(wcs);
829
830 return cpl_error_get_code();
831}
832/* -------------------------------------------------------------------------- */
833
834
835/* -------------------------------------------------------------------------- */
862/*----------------------------------------------------------------------------*/
863/* unit tests utility */
864static cpl_propertylist*
865hdrl_resample_util_crea_header_image(const int naxis, const int sx, const int sy,
866 const double ra, const double dec,
867 const double cd11, const double cd12,
868 const double cd21, const double cd22,
869 const double crpix1, const double crpix2,
870 const double crval1, const double crval2,
871 const double cdelt1, const double cdelt2,
872 const char* ctype1, const char* ctype2,
873 const char* cunit1, const char* cunit2)
874{
875 cpl_propertylist* plist = cpl_propertylist_new();
876 cpl_propertylist_append_int(plist, "NAXIS", naxis);
877
878 cpl_propertylist_append_int(plist, "NAXIS1", sx);
879 cpl_propertylist_append_int(plist, "NAXIS2", sy);
880
881 cpl_propertylist_append_double(plist, "RA", ra);
882 cpl_propertylist_append_double(plist, "DEC", dec);
883
884 cpl_propertylist_append_double(plist, "CRPIX1", crpix1);
885 cpl_propertylist_append_double(plist, "CRPIX2", crpix2);
886
887 cpl_propertylist_append_double(plist, "CRVAL1", crval1);
888 cpl_propertylist_append_double(plist, "CRVAL2", crval2);
889
890 cpl_propertylist_append_double(plist, "CDELT1", cdelt1);
891 cpl_propertylist_append_double(plist, "CDELT2", cdelt2);
892
893 cpl_propertylist_append_string(plist, "CTYPE1", ctype1);
894 cpl_propertylist_append_string(plist, "CTYPE2", ctype2);
895
896 cpl_propertylist_append_string(plist, "CUNIT1", cunit1);
897 cpl_propertylist_append_string(plist, "CUNIT2", cunit2);
898
899 cpl_propertylist_append_double(plist, "CD1_1", cd11);
900 cpl_propertylist_append_double(plist, "CD1_2", cd12);
901 cpl_propertylist_append_double(plist, "CD2_1", cd21);
902 cpl_propertylist_append_double(plist, "CD2_2", cd22);
903
904 /* To be sure to have a standard FITS header we save and reload the image */
905 cpl_image* ima = cpl_image_new(sx, sy, CPL_TYPE_INT);
906 cpl_image_add_scalar(ima, 1);
907
908 cpl_image_save(ima, "ima.fits", CPL_TYPE_INT, plist, CPL_IO_DEFAULT);
909 cpl_image_delete(ima);
910 cpl_propertylist_delete(plist);
911 plist = cpl_propertylist_load("ima.fits", 0);
912 cpl_test_error(CPL_ERROR_NONE);
913 return plist;
914}
915
916
961/* unit tests utility */
962static cpl_propertylist*
963hdrl_resample_crea_header_cube(const int naxis,
964 const int sx, const int sy, const int sz,
965 const double ra, const double dec,
966 const double cd11, const double cd12,
967 const double cd21, const double cd22,
968 const double cd13, const double cd31,
969 const double cd23, const double cd32,
970 const double cd33,
971 const double crpix1, const double crpix2, const double crpix3,
972 const double crval1, const double crval2, const double crval3,
973 const double cdelt1, const double cdelt2, const double cdelt3,
974 const char* ctype1, const char* ctype2, const char* ctype3,
975 const char* cunit1, const char* cunit2, const char* cunit3)
976{
977 /* crea first the FITS header for a 2D example */
978 cpl_propertylist* plist =
979 hdrl_resample_util_crea_header_image(naxis, sx, sy, ra, dec,
980 cd11, cd12, cd21, cd22, crpix1, crpix2,
981 crval1, crval2, cdelt1, cdelt2, ctype1,
982 ctype2, cunit1, cunit2);
983
984 /* then add information for a 3D example */
985 cpl_propertylist_update_int(plist, "NAXIS", naxis);
986
987 cpl_propertylist_append_int(plist, "NAXIS3", sz);
988
989 cpl_propertylist_append_double(plist, "CRVAL3", crval3);
990 cpl_propertylist_append_double(plist, "CRPIX3", crpix3);
991 cpl_propertylist_append_double(plist, "CDELT3", cdelt3);
992 cpl_propertylist_append_string(plist, "CTYPE3", ctype3);
993 cpl_propertylist_append_string(plist, "CUNIT3", cunit3);
994
995 cpl_propertylist_append_double(plist, "CD1_3", cd13);
996 cpl_propertylist_append_double(plist, "CD3_1", cd31);
997 cpl_propertylist_append_double(plist, "CD2_3", cd23);
998 cpl_propertylist_append_double(plist, "CD3_2", cd32);
999 cpl_propertylist_append_double(plist, "CD3_3", cd33);
1000
1001
1002 /* To be sure to have a standard FITS header we save & reload the imagelist */
1003 cpl_image* ima = cpl_image_new(sx, sy, CPL_TYPE_INT);
1004 cpl_image_add_scalar(ima, 1);
1005 cpl_imagelist* iml = cpl_imagelist_new();
1006 for(int i = 0; i< sz;i++) {
1007 cpl_imagelist_set(iml, ima, i);
1008 }
1009
1010 cpl_imagelist_save(iml, "iml.fits", CPL_TYPE_INT, plist, CPL_IO_DEFAULT);
1011 cpl_imagelist_delete(iml);
1012 cpl_propertylist_delete(plist);
1013 plist = cpl_propertylist_load("iml.fits", 0);
1014 cpl_test_error(CPL_ERROR_NONE);
1015 return plist;
1016}
1017
1021static cpl_error_code
1022test_invalid_input_hdrl_resample_weight_function_renka(void){
1023 double r;
1024 double r_c;
1025
1026 /* test invalid input */
1027 /* case r = 0 */
1028 r = 0;
1029 r_c = 1;
1030 cpl_test_abs(hdrl_resample_weight_function_renka(r, r_c), FLT_MAX, HDRL_EPS_TEST);
1031 cpl_test_error(CPL_ERROR_NONE);
1032
1033 /* case r > r_c */
1034 r = 2;
1035 r_c = 1;
1036 cpl_test_abs(hdrl_resample_weight_function_renka(r, r_c), DBL_MIN, HDRL_EPS_TEST);
1037 cpl_test_error(CPL_ERROR_NONE);
1038
1039 /* case r == r_c */
1040 r = 2;
1041 r_c = 2;
1042 cpl_test_abs(hdrl_resample_weight_function_renka(r, r_c), DBL_MIN, HDRL_EPS_TEST);
1043 cpl_test_error(CPL_ERROR_NONE);
1044 return cpl_error_get_code();
1045}
1046static cpl_error_code
1047test_hdrl_resample_weight_function_renka(void)
1048{
1049 double r;
1050 double r_c;
1051 double res;
1052
1053 /* test invalid input */
1054 test_invalid_input_hdrl_resample_weight_function_renka();
1055
1056 /* test valid input */
1057 /* case r < r_c */
1058 r = 2.;
1059 r_c = 3.;
1060
1061 res = 0.0277777777777777762; /* value from print out on fc32 (laptop)*/
1062 cpl_test_abs(hdrl_resample_weight_function_renka(r, r_c), res, HDRL_EPS_TEST);
1063
1064 cpl_test_error(CPL_ERROR_NONE);
1065
1066 return cpl_error_get_code();
1067}
1068
1072
1073static cpl_error_code
1074test_hdrl_resample_weight_function_drizzle(void)
1075{
1076 double x_in, y_in, z_in;
1077 double x_out, y_out, z_out;
1078 double dx, dy, dz;
1079
1080 /* 1st test case */
1081 x_in = 2, y_in = 2, z_in = 2;
1082 dx = 1, dy = 1, dz = 1;
1083 x_out = 1, y_out = 1, z_out = 1;
1084 double res = 0.015625;/* value from computation in code */
1085 hdrl_resample_weight_function_drizzle(x_in, y_in, z_in, x_out, y_out, z_out,
1086 dx, dy, dz);
1087 cpl_test_error(CPL_ERROR_NONE);
1088
1089 cpl_test_abs(hdrl_resample_weight_function_drizzle(x_in, y_in, z_in,
1090 x_out, y_out, z_out, dx, dy, dz),res,HDRL_EPS_TEST);
1091
1092 /* 2nd test case */
1093 x_in = 3, y_in = 3, z_in = 3;
1094 dx = 1, dy = 1, dz = 1;
1095 x_out = 2, y_out = 2, z_out = 2;
1096 res = 0.125;/* value from computation in code */
1097 hdrl_resample_weight_function_drizzle(x_in, y_in, z_in, x_out, y_out, z_out,
1098 dx, dy, dz);
1099 cpl_test_error(CPL_ERROR_NONE);
1100 cpl_test_abs(hdrl_resample_weight_function_drizzle(x_in, y_in, z_in,
1101 x_out, y_out, z_out, dx, dy, dz),res,HDRL_EPS_TEST);
1102
1103 /* 3rd test case */
1104 x_in = 2, y_in = 2, z_in = 2;
1105 dx = 10, dy = 10, dz = 10;
1106 x_out = 1, y_out = 1, z_out = 1;
1107 res = 0; /* value from computation in code */
1108 cpl_test_abs(hdrl_resample_weight_function_drizzle(x_in, y_in, z_in,
1109 x_out, y_out, z_out, dx, dy, dz),res,HDRL_EPS_TEST);
1110
1111 cpl_test_error(CPL_ERROR_NONE);
1112
1113 return cpl_error_get_code();
1114}
1115
1119static cpl_error_code
1120test_hdrl_resample_weight_function_linear(void)
1121{
1122 /* test invalid input */
1123 /* r = 0 */
1124 cpl_test_abs(hdrl_resample_weight_function_linear(0), FLT_MAX, HDRL_EPS_TEST);
1125 cpl_test_error(CPL_ERROR_NONE);
1126
1127 /* test valid input */
1128 /* r = 2 */
1129 cpl_test_abs(hdrl_resample_weight_function_linear(2), 0.5, HDRL_EPS_TEST);
1130 cpl_test_error(CPL_ERROR_NONE);
1131
1132
1133 return cpl_error_get_code();
1134}
1135
1139static cpl_error_code
1140test_hdrl_resample_weight_function_quadratic(void)
1141{
1142 /* test invalid input */
1143 /* r = 0 */
1144 cpl_test_abs(hdrl_resample_weight_function_quadratic(0), FLT_MAX, HDRL_EPS_TEST);
1145 cpl_test_error(CPL_ERROR_NONE);
1146
1147 /* test valid input */
1148 /* r = 4 */
1149 cpl_test_abs(hdrl_resample_weight_function_quadratic(4), 0.25, HDRL_EPS_TEST);
1150 cpl_test_error(CPL_ERROR_NONE);
1151
1152
1153 return cpl_error_get_code();
1154}
1155
1159static cpl_error_code
1160test_hdrl_resample_weight_function_sinc(void)
1161{
1162
1163 double r;
1164 double res;
1165 /* test invalid input */
1166 /* r = 0 */
1167 cpl_test_abs(hdrl_resample_weight_function_sinc(DBL_EPSILON), 1., HDRL_EPS_TEST);
1168 cpl_test_error(CPL_ERROR_NONE);
1169
1170 /* test valid input */
1171 r = 0.25;
1172
1173 res = 0.900316316157106056;/* value from computation in code */
1174 cpl_test_abs(hdrl_resample_weight_function_sinc(r), res, HDRL_EPS_TEST);
1175 cpl_test_error(CPL_ERROR_NONE);
1176
1177 return cpl_error_get_code();
1178}
1179
1183static cpl_error_code
1184test_hdrl_resample_weight_function_lanczos(void)
1185{
1186
1187 double dx = 1.;
1188 double dy = 1.;
1189 double dz = 1.;
1190 unsigned int ld = 4; /* loop distance */
1191 unsigned int lks = 2; /* lanczos kernel size */
1192
1193 /* test invalid input */
1194 /* dx = 0 || dy =0 || dz = 0 */
1195 cpl_test_abs(hdrl_resample_weight_function_lanczos(5., dy, dz, ld, lks), 0., HDRL_EPS_TEST);
1196 cpl_test_abs(hdrl_resample_weight_function_lanczos(dx, 5., dz, ld, lks), 0., HDRL_EPS_TEST);
1197 cpl_test_abs(hdrl_resample_weight_function_lanczos(dx, dy, 5., ld, lks), 0., HDRL_EPS_TEST);
1198 cpl_test_error(CPL_ERROR_NONE);
1199
1200 /* test valid input */
1201 double res = 4.32283142061004719e-50;/* value from computation in code */
1202 cpl_test_abs(hdrl_resample_weight_function_lanczos(dx, dy, dz, ld, lks), res, HDRL_EPS_TEST);
1203 cpl_test_error(CPL_ERROR_NONE);
1204
1205 return cpl_error_get_code();
1206}
1207
1211static cpl_image*
1212hdrl_resample_util_crea_image_dice_5(const cpl_size sx, const cpl_size sy,
1213 const double value)
1214{
1215 assert(value > 0);
1216 cpl_image * image = cpl_image_new(sx, sy, CPL_TYPE_DOUBLE);
1217 cpl_image * dice = cpl_image_new(sx, sy, CPL_TYPE_DOUBLE);
1218 cpl_image_add_scalar(image,value);
1219
1220 cpl_size xc = 0.5 * sx;
1221 cpl_size yc = 0.5 * sy;
1222 cpl_size dx = 2;
1223 cpl_size dy = 2;
1224 cpl_size xl = 0.25 * sx;
1225 cpl_size yl = 0.25 * sy;
1226 cpl_size xh = 0.75 * sx;
1227 cpl_size yh = 0.75 * sy;
1228 double flux_point=2.*value;
1229 cpl_image_fill_window(dice, xc - dx, yc - dy, xc + dx, yc + dy, flux_point);
1230 cpl_image_fill_window(dice, xl - dx, yl - dy, xl + dx, yl + dy, flux_point);
1231 cpl_image_fill_window(dice, xl - dx, yh - dy, xl + dx, yh + dy, flux_point);
1232 cpl_image_fill_window(dice, xh - dx, yh - dy, xh + dx, yh + dy, flux_point);
1233 cpl_image_fill_window(dice, xh - dx, yl - dy, xh + dx, yl + dy, flux_point);
1234 cpl_image_add(image,dice);
1235
1236 cpl_image_delete(dice);
1237 cpl_test_error(CPL_ERROR_NONE);
1238 return image;
1239}
1240
1241
1242/*----------------------------------------------------------------------------*/
1254/*----------------------------------------------------------------------------*/
1255static hdrl_parameter *
1256hdrl_resample_util_methodparam_create (
1257 const int loop_distance,
1258 const double critical_radius_renka,
1259 const int kernel_size_lanczos,
1260 const double pix_frac_drizzle_x,
1261 const double pix_frac_drizzle_y,
1262 const double pix_frac_drizzle_l,
1263 const int resample_method,
1264 const cpl_boolean use_errorweights
1265)
1266{
1267 hdrl_parameter *aParams_method = NULL;
1268
1269 /* Create the right re-sampling parameter */
1270 switch(resample_method) {
1271 case HDRL_RESAMPLE_METHOD_NEAREST:
1272 aParams_method =
1273 hdrl_resample_parameter_create_nearest();
1274 break;
1275 case HDRL_RESAMPLE_METHOD_RENKA:
1276 aParams_method =
1277 hdrl_resample_parameter_create_renka(loop_distance, use_errorweights,
1278 critical_radius_renka);
1279 break;
1280 case HDRL_RESAMPLE_METHOD_LINEAR:
1281 aParams_method =
1282 hdrl_resample_parameter_create_linear(loop_distance, use_errorweights);
1283 break;
1284 case HDRL_RESAMPLE_METHOD_QUADRATIC:
1285 aParams_method =
1286 hdrl_resample_parameter_create_quadratic(loop_distance, use_errorweights);
1287 break;
1288 case HDRL_RESAMPLE_METHOD_DRIZZLE:
1289 aParams_method =
1290 hdrl_resample_parameter_create_drizzle(loop_distance, use_errorweights,
1291 pix_frac_drizzle_x,
1292 pix_frac_drizzle_y,
1293 pix_frac_drizzle_l);
1294 break;
1295 case HDRL_RESAMPLE_METHOD_LANCZOS:
1296 aParams_method =
1297 hdrl_resample_parameter_create_lanczos(loop_distance, use_errorweights,
1298 kernel_size_lanczos);
1299 break;
1300 default:
1301 aParams_method =
1302 hdrl_resample_parameter_create_lanczos(loop_distance, use_errorweights,
1303 kernel_size_lanczos);
1304 cpl_msg_warning (cpl_func,
1305 "You set an unsupported method! Default to LANCZOS");
1306 break;
1307 }
1308
1309 cpl_test_error(CPL_ERROR_NONE);
1310 return aParams_method;
1311}
1312
1313
1314/*----------------------------------------------------------------------------*/
1321/*----------------------------------------------------------------------------*/
1322static cpl_wcs *
1323hdrl_resample_util_wcs_create_example_params(void)
1324{
1325
1326 cpl_propertylist* p = cpl_propertylist_new();
1327 /* leave out the velocity and stokes axes */
1328 cpl_propertylist_append_int(p, "NAXIS", 2);
1329 cpl_propertylist_append_int(p, "NAXIS1", 512);
1330 cpl_propertylist_append_int(p, "NAXIS2", 512);
1331 cpl_propertylist_append_double(p, "CRPIX1", 256.);
1332
1333 /* use the CDi_j matrix instead of CDELT */
1334 cpl_propertylist_append_double(p, "CD1_1", -0.003);
1335 cpl_propertylist_append_string(p, "CTYPE1", "RA---TAN");
1336 cpl_propertylist_append_double(p, "CRVAL1", 45.83);
1337 cpl_propertylist_append_string(p, "CUNIT1", "deg");
1338
1339 cpl_propertylist_append_double(p, "CRPIX2", 257.);
1340 cpl_propertylist_append_double(p, "CD2_2", 0.003);
1341 cpl_propertylist_append_string(p, "CTYPE2", "DEC--TAN");
1342 cpl_propertylist_append_double(p, "CRVAL2", 63.57);
1343 cpl_propertylist_append_string(p, "CUNIT2", "deg");
1344
1345 /* no cross terms */
1346 cpl_propertylist_append_double(p, "CD1_2", 0.);
1347 cpl_propertylist_append_double(p, "CD2_1", 0.);
1348
1349 /* XXX RA,DEC identical to the CRVALi for our functions */
1350 cpl_propertylist_append_double(p, "RA", 45.83);
1351 cpl_propertylist_append_double(p, "DEC", 63.57);
1352
1353 cpl_wcs* wcs = cpl_wcs_new_from_propertylist(p);
1354 cpl_propertylist_delete(p);
1355 cpl_test_error(CPL_ERROR_NONE);
1356 return wcs;
1357} /* hdrl_resample_util_wcs_create_example_params() */
1358
1359/*----------------------------------------------------------------------------*/
1364/*----------------------------------------------------------------------------*/
1365static cpl_wcs*
1366hdrl_resample_util_crea_wcs_2d(void)
1367{
1368 cpl_propertylist* plist = NULL;
1369 /* test valid input */
1370 int sx = HDRL_SIZE_X;
1371 int sy = HDRL_SIZE_Y;
1372
1373 int naxis = 2;
1374 double ra = 10.;
1375 double dec = 20.;
1376 double cd11 = HDRL_CD11;
1377 double cd22 = HDRL_CD22;
1378 double cd12 = HDRL_CD12;
1379 double cd21 = HDRL_CD21;
1380 const double crpix1 = HDRL_CRPIX1;
1381 const double crpix2 = HDRL_CRPIX2;
1382 const double crval1 = HDRL_CRVAL1;
1383 const double crval2 = HDRL_CRVAL2;
1384 const double cdelt1 = HDRL_CDELT1;
1385 const double cdelt2 = HDRL_CDELT2;
1386 const char* cunit1 = "";
1387 const char* cunit2 = "";
1388 const char* ctype1 = "pix";
1389 const char* ctype2 = "pix";
1390
1391 plist = hdrl_resample_util_crea_header_image(naxis, sx, sy, ra, dec,
1392 cd11, cd12, cd21, cd22, crpix1,
1393 crpix2, crval1, crval2,
1394 cdelt1, cdelt2, ctype1, ctype2,
1395 cunit1, cunit2);
1396 cpl_test_error(CPL_ERROR_NONE);
1397
1398 cpl_wcs* wcs = cpl_wcs_new_from_propertylist(plist);
1399 cpl_test_nonnull(wcs);
1400 cpl_test_error(CPL_ERROR_NONE);
1401 cpl_propertylist_delete(plist);
1402 return wcs;
1403
1404}
1405/*----------------------------------------------------------------------------*/
1410/*----------------------------------------------------------------------------*/
1411static cpl_wcs*
1412hdrl_resample_util_crea_wcs_3d(void)
1413{
1414 /* test 3D case */
1415 cpl_propertylist* plist = NULL;
1416 int sx = HDRL_SIZE_X;
1417 int sy = HDRL_SIZE_Y;
1418
1419 int naxis = 2;
1420 double ra = 10.;
1421 double dec = 20.;
1422 double cd11 = HDRL_CD11;
1423 double cd22 = HDRL_CD22;
1424 double cd12 = HDRL_CD12;
1425 double cd21 = HDRL_CD21;
1426 const double crpix1 = HDRL_CRPIX1;
1427 const double crpix2 = HDRL_CRPIX2;
1428 const double crval1 = HDRL_CRVAL1;
1429 const double crval2 = HDRL_CRVAL2;
1430 const double cdelt1 = HDRL_CDELT1;
1431 const double cdelt2 = HDRL_CDELT2;
1432 const char* cunit1 = "deg";
1433 const char* cunit2 = "deg";
1434 const char* ctype1 = "RA---TAN";
1435 const char* ctype2 = "DEC--TAN";
1436
1437 int sz = HDRL_SIZE_Z;
1438 double cd13 = HDRL_CD13;
1439 double cd31 = HDRL_CD31;
1440 double cd23 = HDRL_CD23;
1441 double cd32 = HDRL_CD32;
1442 double cd33 = HDRL_CD33;
1443 const double crpix3 = HDRL_CRPIX3;
1444 const double crval3 = HDRL_CRVAL3;
1445 const double cdelt3 = HDRL_CDELT3;
1446 const char* cunit3 = "m";
1447 const char* ctype3 = "WAV";
1448 naxis = 3;
1449 cpl_test_error(CPL_ERROR_NONE);
1450 cpl_propertylist_delete(plist);
1451 cpl_test_error(CPL_ERROR_NONE);
1452 plist = hdrl_resample_crea_header_cube(naxis, sx, sy, sz, ra, dec,
1453 cd11, cd12, cd21, cd22, cd13, cd31,
1454 cd23, cd32, cd33,
1455 crpix1, crpix2, crpix3,
1456 crval1, crval2, crval3,
1457 cdelt1, cdelt2, cdelt3,
1458 ctype1, ctype2, ctype3,
1459 cunit1, cunit2, cunit3);
1460 cpl_test_error(CPL_ERROR_NONE);
1461
1462 cpl_wcs* wcs = cpl_wcs_new_from_propertylist(plist);
1463 cpl_test_nonnull(wcs);
1464 cpl_test_error(CPL_ERROR_NONE);
1465 cpl_propertylist_delete(plist);
1466 return wcs;
1467
1468}
1469
1470CPL_DIAG_PRAGMA_PUSH_IGN(-Wunused-function)
1471/*----------------------------------------------------------------------------*/
1476/*----------------------------------------------------------------------------*/
1477static cpl_error_code
1478test_resample_wcs_as_muse(void)
1479{
1480
1481 // const double kLimitPos = 0.05; /* [pix] object position detection accuracy */
1482 /* accuracy limits for coordinate transformation tests */
1483 // const double kLimitDeg = DBL_EPSILON * 115; /* ~10 nano-arcsec for trafo to deg */
1484 // const double kLimitDegF = FLT_EPSILON * 13.51; /* ~5.8 milli-arcsec for trafo to *
1485 // * deg, with value stored in float */
1486 // const double kLimitPix = FLT_EPSILON * 205; /* ~1/40000th pixel for trafo to pix */
1487 // const double kLimitRot = FLT_EPSILON * 10; /* ~4.3 milli-arcsec for rotations */
1488 // const double kLimitSca = DBL_EPSILON * 7461; /* pixel scales to ~5.5 nano-arcsec */
1489
1490 const double kLimitPPl = FLT_EPSILON * 2.88; /* ~1.24 milli-arcsec in proj. plane */
1491
1492
1493 double v1, v2;
1494 /***************************************************************************
1495 * test transformations between projection plane and celestial coordinates *
1496 ***************************************************************************/
1497 /* use values from WCS Paper II, example 1 as references */
1498
1499 cpl_test_error(CPL_ERROR_NONE);
1500 cpl_wcs* wcs = hdrl_resample_util_wcs_create_example_params();
1501 cpl_test_nonnull(wcs);
1502 /* NB: for the following we just want to test the WCS==> other coordinates
1503 * are casual.
1504 */
1505 hdrl_parameter* hpar =
1506 hdrl_resample_parameter_create_outgrid2D_userdef(1, 1, 10, 10.1, 10,
1507 10.1, 5.);
1508 cpl_test_nonnull(hpar);
1509 hdrl_resample_outgrid_parameter * par = (hdrl_resample_outgrid_parameter *) hpar;
1510 cpl_test_nonnull(par);
1511 par->wcs = wcs;
1512 double x, y;
1513 cpl_errorstate state = cpl_errorstate_get();
1514 cpl_test_error(CPL_ERROR_NONE);
1515
1516 cpl_test(hdrl_resample_wcs_projplane_from_celestial(par, 47.503264, 62.795111,
1517 &x, &y) == CPL_ERROR_NONE);
1518 cpl_test_error(CPL_ERROR_NONE);
1519 cpl_test(fabs(x - 0.765000) < kLimitPPl && fabs(y - (-0.765000)) < kLimitPPl);
1520 cpl_msg_debug(__func__, "SE corner: %f,%f (%e,%e <? %e", x, y, x - 0.765000,
1521 y - (-0.765000), kLimitPPl);
1522 cpl_test(hdrl_resample_wcs_projplane_from_celestial(par, 47.595581, 64.324332,
1523 &x, &y) == CPL_ERROR_NONE);
1524 cpl_test(fabs(x - 0.765000) < kLimitPPl && fabs(y - 0.765000) < kLimitPPl);
1525 cpl_msg_debug(__func__, "NE corner: %f,%f (%e,%e <? %e", x, y, x - 0.765000,
1526 y - 0.765000, kLimitPPl);
1527 cpl_test(hdrl_resample_wcs_projplane_from_celestial(par, 44.064419, 64.324332,
1528 &x, &y) == CPL_ERROR_NONE);
1529 cpl_test(fabs(x - (-0.765000)) < kLimitPPl && fabs(y - 0.765000) < kLimitPPl);
1530 cpl_msg_debug(__func__, "NW corner: %f,%f (%e,%e <? %e", x, y, x - (-0.765000),
1531 y - 0.765000, kLimitPPl);
1532 state = cpl_errorstate_get();
1533 cpl_test(hdrl_resample_wcs_projplane_from_celestial(NULL, 1., 1., &v1, &v2)
1534 == CPL_ERROR_NULL_INPUT);
1535 cpl_test(hdrl_resample_wcs_projplane_from_celestial(par, 1., 1., NULL, &v2)
1536 == CPL_ERROR_NULL_INPUT);
1537 cpl_test(hdrl_resample_wcs_projplane_from_celestial(par, 1., 1., &v1, NULL)
1538 == CPL_ERROR_NULL_INPUT);
1539 cpl_errorstate_set(state);
1540
1541 cpl_test_error(CPL_ERROR_NONE);
1542 cpl_test(hdrl_resample_wcs_projplane_from_celestial(par, 1., 1., &v1, &v2)
1543 == CPL_ERROR_NONE);
1544 cpl_test_error(CPL_ERROR_NONE);
1545 cpl_errorstate_set(state);
1546
1547 cpl_wcs_delete(wcs);
1548 hdrl_parameter_destroy(hpar);
1549 cpl_test_error(CPL_ERROR_NONE);
1550 return cpl_error_get_code();
1551}
1552
1553CPL_DIAG_PRAGMA_POP
1554/*----------------------------------------------------------------------------*/
1559/*----------------------------------------------------------------------------*/
1560static cpl_error_code
1561test_hdrl_resample_image_to_table(void)
1562{
1563
1564 cpl_table* tab = NULL;
1565 cpl_image* data = NULL;
1566 cpl_image* error = NULL;
1567 cpl_image* quality = NULL;
1568 cpl_wcs* wcs = NULL;
1569
1570 /* test invalid input */
1571 tab = hdrl_resample_image_to_table(NULL, wcs);
1572 cpl_test_error(CPL_ERROR_NULL_INPUT);
1573 cpl_test_null(tab);
1574
1575 /* test valid input */
1576 int sx = HDRL_SIZE_X;
1577 int sy = HDRL_SIZE_Y;
1578 double value = HDRL_FLUX_ADU;
1579
1580 data = cpl_image_new(sx, sy, CPL_TYPE_DOUBLE);
1581 cpl_image_add_scalar(data, value);
1582
1583 error = cpl_image_power_create(data, 0.5);
1584 quality = cpl_image_new(sx, sy, CPL_TYPE_INT);
1585 int naxis = 2;
1586 double ra = 10.;
1587 double dec = 20.;
1588 const double cd11 = -3.47222e-05;
1589 const double cd22 = 3.47222e-05;
1590 const double cd12 = 0;
1591 const double cd21 = 0;
1592 const double crpix1 = 33.5;
1593 const double crpix2 = 33.5;
1594 const double crval1 = 48.0718057375143246;
1595 const double crval2 = -20.6230284673176705;
1596 const double cdelt1 = 0;
1597 const double cdelt2 = 0;
1598 const char* cunit1 = "deg";
1599 const char* cunit2 = "deg";
1600 const char* ctype1 = "RA---TAN";
1601 const char* ctype2 = "DEC--TAN";
1602
1603 cpl_propertylist* plist =
1604 hdrl_resample_util_crea_header_image(naxis, sx, sy, ra, dec,
1605 cd11, cd12, cd21, cd22,
1606 crpix1, crpix2, crval1, crval2,
1607 cdelt1, cdelt2, ctype1, ctype2,
1608 cunit1, cunit2);
1609
1610 wcs = cpl_wcs_new_from_propertylist(plist);
1611 cpl_test_nonnull(wcs);
1612 hdrl_image* hima = hdrl_resample_util_hdrl_image_create(data, error, quality);
1613 tab = hdrl_resample_image_to_table(hima, wcs);
1614 cpl_test_error(CPL_ERROR_NONE);
1615
1616 /* test expected table columns */
1617 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_RA));
1618 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_DEC));
1619 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_LAMBDA));
1620 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_BPM));
1621 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_ERRORS));
1622
1623 double* ptablambda = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_LAMBDA);
1624 double* ptabdata = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_DATA);
1625 int* ptabbpm = cpl_table_get_data_int(tab, HDRL_RESAMPLE_TABLE_BPM);
1626 double* ptaberr = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_ERRORS);
1627
1628 cpl_test_abs(0, ptablambda[0], HDRL_EPS_TEST);
1629 cpl_test_abs(HDRL_FLUX_ADU, ptabdata[0], HDRL_EPS_TEST);
1630 cpl_test_abs(0, ptabbpm[0], HDRL_EPS_TEST);
1631 cpl_test_abs(10, ptaberr[0], HDRL_EPS_TEST);
1632
1633 /* free memory */
1634 cpl_table_delete(tab);
1635 cpl_image_delete(data);
1636 cpl_image_delete(error);
1637 cpl_image_delete(quality);
1638 hdrl_image_delete(hima);
1639 cpl_propertylist_delete(plist);
1640 cpl_wcs_delete(wcs);
1641 cpl_test_error(CPL_ERROR_NONE);
1642 return cpl_error_get_code();
1643}
1644/*----------------------------------------------------------------------------*/
1649/*----------------------------------------------------------------------------*/
1650static cpl_error_code
1651test_hdrl_resample_imagelist_to_table(void)
1652{
1653 /* test invalid input */
1654 cpl_table* tab = NULL;
1655 cpl_imagelist* ilist = NULL;
1656 cpl_imagelist* elist = NULL;
1657 cpl_imagelist* qlist = NULL;
1658 cpl_image* data = NULL;
1659 cpl_image* error = NULL;
1660 cpl_image* quality = NULL;
1661 cpl_wcs* wcs = NULL;
1662
1663 /* test invalid input */
1664 tab = hdrl_resample_imagelist_to_table(NULL, wcs);
1665 cpl_test_error(CPL_ERROR_NULL_INPUT);
1666 cpl_test_null(tab);
1667
1668 /* test valid input */
1669 int sx = 67;
1670 int sy = 67;
1671 double value = HDRL_FLUX_ADU;
1672
1673 data = cpl_image_new(sx, sy, CPL_TYPE_DOUBLE);
1674 cpl_image_add_scalar(data, value);
1675
1676 error = cpl_image_power_create(data, 0.5);
1677 quality = cpl_image_new(sx, sy, CPL_TYPE_INT);
1678 int naxis = 3;
1679 double ra = 10.;
1680 double dec = 20.;
1681 const double cd11 = -3.47222e-05;
1682 const double cd22 = 3.47222e-05;
1683 const double cd12 = 0;
1684 const double cd21 = 0;
1685 const double crpix1 = 33.5;
1686 const double crpix2 = 33.5;
1687 const double crval1 = 48.0706;
1688 const double crval2 = -20.6219;
1689 const double cdelt1 = 0;
1690 const double cdelt2 = 0;
1691 const char* cunit1 = "deg";
1692 const char* cunit2 = "deg";
1693 const char* ctype1 = "RA---TAN";
1694 const char* ctype2 = "DEC--TAN";
1695 int sz = 2218;
1696 double cd13 = 0;
1697 double cd31 = 0;
1698 double cd23 = 0;
1699 double cd32 = 0;
1700 double cd33 = 2.45e-10;
1701 const double crpix3 = 1;
1702 const double crval3 = 1.9283e-06;
1703 const double cdelt3 = 0.1;
1704 const char* cunit3 = "m";
1705 const char* ctype3 = "WAVE";
1706
1707 cpl_test_error(CPL_ERROR_NONE);
1708 cpl_propertylist* plist =
1709 hdrl_resample_crea_header_cube(naxis, sx, sy, sz, ra, dec,
1710 cd11, cd12, cd21, cd22,
1711 cd13, cd31, cd23, cd32, cd33,
1712 crpix1, crpix2, crpix3,
1713 crval1, crval2, crval3,
1714 cdelt1, cdelt2, cdelt3,
1715 ctype1, ctype2, ctype3,
1716 cunit1, cunit2, cunit3);
1717
1718 ilist = cpl_imagelist_new();
1719 elist = cpl_imagelist_new();
1720 qlist = cpl_imagelist_new();
1721
1722 cpl_imagelist_set(ilist, data, 0);
1723 cpl_imagelist_set(elist, error, 0);
1724 cpl_imagelist_set(qlist, quality, 0);
1725
1726 wcs = cpl_wcs_new_from_propertylist(plist);
1727 cpl_test_nonnull(wcs);
1728 cpl_test_error(CPL_ERROR_NONE);
1729 hdrl_imagelist* hlist = hdrl_resample_util_hdrl_imagelist_create(ilist, elist,
1730 qlist);
1731 tab = hdrl_resample_imagelist_to_table(hlist, wcs);
1732 cpl_test_error(CPL_ERROR_NONE);
1733
1734 /* test expected table columns */
1735 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_RA));
1736 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_DEC));
1737 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_LAMBDA));
1738 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_BPM));
1739 cpl_test_eq(1, cpl_table_has_column(tab, HDRL_RESAMPLE_TABLE_ERRORS));
1740
1741 //double* ptabxpos = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_RA);
1742 //double* ptabypos = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_DEC);
1743 double* ptablambda = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_LAMBDA);
1744 double* ptabdata = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_DATA);
1745 int* ptabbpm = cpl_table_get_data_int(tab, HDRL_RESAMPLE_TABLE_BPM);
1746 double* ptaberr = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_ERRORS);
1747
1748 /* test expected values */
1749 cpl_test_abs(0, ptablambda[0], HDRL_EPS_TEST);
1750 cpl_test_abs(HDRL_FLUX_ADU, ptabdata[0], HDRL_EPS_TEST);
1751 cpl_test_abs(0, ptabbpm[0], HDRL_EPS_TEST);
1752 cpl_test_abs(10, ptaberr[0], HDRL_EPS_TEST);
1753
1754 /* free memory */
1755 cpl_imagelist_delete(ilist);
1756 cpl_imagelist_delete(elist);
1757 cpl_imagelist_delete(qlist);
1758 hdrl_imagelist_delete(hlist);
1759 cpl_table_delete(tab);
1760 cpl_propertylist_delete(plist);
1761 cpl_wcs_delete(wcs);
1762 cpl_test_error(CPL_ERROR_NONE);
1763 return cpl_error_get_code();
1764}
1765/*----------------------------------------------------------------------------*/
1770/*----------------------------------------------------------------------------*/
1771static cpl_error_code
1772test_hdrl_wcs_to_propertylist(void) {
1773
1774 cpl_propertylist* plist = NULL;
1775 cpl_wcs* wcs = NULL;
1776 cpl_boolean only2d = CPL_TRUE;
1777
1778 /* test improper input */
1779 hdrl_wcs_to_propertylist(NULL, plist, only2d);
1780 cpl_test_error(CPL_ERROR_NULL_INPUT);
1781
1782 hdrl_wcs_to_propertylist(wcs, NULL, only2d);
1783 cpl_test_error(CPL_ERROR_NULL_INPUT);
1784
1785 /* test proper input */
1786 wcs = hdrl_resample_util_crea_wcs_2d();
1787 cpl_test_error(CPL_ERROR_NONE);
1788 plist = cpl_propertylist_new();
1789 hdrl_wcs_to_propertylist(wcs, plist, only2d);
1790 cpl_test_error(CPL_ERROR_NONE);
1791
1792 cpl_wcs_delete(wcs);
1793 wcs = hdrl_resample_util_crea_wcs_3d();
1794 cpl_test_error(CPL_ERROR_NONE);
1795
1796 hdrl_wcs_to_propertylist(wcs, plist, !only2d);
1797 cpl_test_error(CPL_ERROR_NONE);
1798
1799 /* start detailed checks on expected values */
1800 const cpl_array *crval = cpl_wcs_get_crval(wcs);
1801 const cpl_array *crpix = cpl_wcs_get_crpix(wcs);
1802 const cpl_array *ctype = cpl_wcs_get_ctype(wcs);
1803 const cpl_array *cunit = cpl_wcs_get_cunit(wcs);
1804 const cpl_matrix *cd = cpl_wcs_get_cd(wcs);
1805 int naxis_out = cpl_wcs_get_image_naxis(wcs);
1806
1807 int naxis = 3;
1808 int sx = HDRL_SIZE_X;
1809 int sz = HDRL_SIZE_Z;
1810 cpl_test_eq(naxis, naxis_out);
1811
1812 /* Check NAXIS */
1813 for (int i = 0; i < naxis_out; i++) {
1814
1815 char * buf = cpl_sprintf("NAXIS%d", i + 1);
1816 /* sx = sy = HDRL_SIZE_X */
1817 if(i<2) {
1818 cpl_test_eq(sx, cpl_propertylist_get_int(plist, buf));
1819 } else {
1820 cpl_test_eq(sz, cpl_propertylist_get_int(plist, buf));
1821 }
1822 cpl_free(buf);
1823 }
1824
1825 int err = 0;
1826 /* CRVAL */
1827 cpl_test_abs(cpl_propertylist_get_double(plist, "CRVAL1"),
1828 cpl_array_get_double(crval, 0, &err), HDRL_EPS_TEST);
1829
1830 cpl_test_abs(cpl_propertylist_get_double(plist, "CRVAL2"),
1831 cpl_array_get_double(crval, 1, &err), HDRL_EPS_TEST);
1832
1833 cpl_test_abs(cpl_propertylist_get_double(plist, "CRVAL3"),
1834 cpl_array_get_double(crval, 2, &err), HDRL_EPS_TEST);
1835
1836 /* CRPIX */
1837 cpl_test_abs(cpl_propertylist_get_double(plist, "CRPIX1"),
1838 cpl_array_get_double(crpix, 0, &err), HDRL_EPS_TEST);
1839
1840 cpl_test_abs(cpl_propertylist_get_double(plist, "CRPIX2"),
1841 cpl_array_get_double(crpix, 1, &err), HDRL_EPS_TEST);
1842
1843 cpl_test_abs(cpl_propertylist_get_double(plist, "CRPIX3"),
1844 cpl_array_get_double(crpix, 2, &err), HDRL_EPS_TEST);
1845
1846 /* CTYPE */
1847 const char* string_val = NULL;
1848 const char* string_chk = NULL;
1849
1850 string_val = cpl_propertylist_get_string(plist, "CTYPE1");
1851 string_chk = cpl_array_get_string(ctype, 0);
1852 cpl_test_eq(0, strcmp(string_val,string_chk));
1853
1854 string_val = cpl_propertylist_get_string(plist, "CTYPE2");
1855 string_chk = cpl_array_get_string(ctype, 1);
1856 cpl_test_eq(0, strcmp(string_val,string_chk));
1857
1858 string_val = cpl_propertylist_get_string(plist, "CTYPE3");
1859 string_chk = cpl_array_get_string(ctype, 2);
1860 cpl_test_eq(0, strcmp(string_val,string_chk));
1861
1862 /* CUNIT */
1863 string_val = cpl_propertylist_get_string(plist, "CUNIT1");
1864 string_chk = cpl_array_get_string(cunit, 0);
1865 cpl_test_eq(0, strcmp(string_val,string_chk));
1866
1867 string_val = cpl_propertylist_get_string(plist, "CUNIT2");
1868 string_chk = cpl_array_get_string(cunit, 1);
1869 cpl_test_eq(0, strcmp(string_val,string_chk));
1870
1871 string_val = cpl_propertylist_get_string(plist, "CUNIT3");
1872 string_chk = cpl_array_get_string(cunit, 2);
1873 cpl_test_eq(0, strcmp(string_val,string_chk));
1874
1875 /* CD */
1876 double cd11 = 1, cd22 = 1, cd33 = 0.1;
1877 double cd12 = 0, cd21 = 0, cd13 = 0, cd31 = 0, cd23 = 0, cd32 = 0;
1878
1879 cd11 = cpl_matrix_get(cd, 0, 0);
1880 cd12 = cpl_matrix_get(cd, 0, 1);
1881 cd21 = cpl_matrix_get(cd, 1, 0);
1882 cd22 = cpl_matrix_get(cd, 1, 1);
1883 cd13 = cpl_matrix_get(cd, 0, 2);
1884 cd31 = cpl_matrix_get(cd, 2, 0);
1885 cd23 = cpl_matrix_get(cd, 1, 2);
1886 cd32 = cpl_matrix_get(cd, 2, 1);
1887 cd33 = cpl_matrix_get(cd, 2, 2);
1888
1889 cpl_test_abs(cd11, cpl_propertylist_get_double(plist, "CD1_1"), HDRL_EPS_TEST);
1890 cpl_test_abs(cd12, cpl_propertylist_get_double(plist, "CD1_2"), HDRL_EPS_TEST);
1891 cpl_test_abs(cd21, cpl_propertylist_get_double(plist, "CD2_1"), HDRL_EPS_TEST);
1892 cpl_test_abs(cd22, cpl_propertylist_get_double(plist, "CD2_2"), HDRL_EPS_TEST);
1893 cpl_test_abs(cd13, cpl_propertylist_get_double(plist, "CD1_3"), HDRL_EPS_TEST);
1894 cpl_test_abs(cd31, cpl_propertylist_get_double(plist, "CD3_1"), HDRL_EPS_TEST);
1895 cpl_test_abs(cd23, cpl_propertylist_get_double(plist, "CD2_3"), HDRL_EPS_TEST);
1896 cpl_test_abs(cd32, cpl_propertylist_get_double(plist, "CD3_2"), HDRL_EPS_TEST);
1897 cpl_test_abs(cd33, cpl_propertylist_get_double(plist, "CD3_3"), HDRL_EPS_TEST);
1898
1899 cpl_propertylist_delete(plist);
1900 cpl_wcs_delete(wcs);
1901 cpl_test_error(CPL_ERROR_NONE);
1902 return cpl_error_get_code();
1903}
1904
1905static cpl_error_code
1906test_resample_invalid_input_outgrid_param(
1907 const double delta_ra,
1908 const double delta_dec,
1909 const double delta_lambda,
1910 const double ra_min,
1911 const double ra_max,
1912 const double dec_min,
1913 const double dec_max,
1914 const double lambda_min,
1915 const double lambda_max,
1916 const double field_margin) {
1917
1918 hdrl_parameter* pErr = NULL;
1919
1920 pErr = hdrl_resample_parameter_create_outgrid2D(0, delta_dec);
1921 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1922 cpl_test_null(pErr);
1923
1924 pErr = hdrl_resample_parameter_create_outgrid2D(delta_ra, 0);
1925 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1926 cpl_test_null(pErr);
1927
1928 pErr = hdrl_resample_parameter_create_outgrid3D(0, delta_dec, delta_lambda);
1929 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1930 cpl_test_null(pErr);
1931
1932 pErr = hdrl_resample_parameter_create_outgrid3D(delta_ra, 0, delta_lambda);
1933 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1934 cpl_test_null(pErr);
1935
1936 pErr = hdrl_resample_parameter_create_outgrid3D(delta_ra, delta_dec, 0);
1937 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1938 cpl_test_null(pErr);
1939
1940 pErr =
1941 hdrl_resample_parameter_create_outgrid2D_userdef(0, delta_dec, ra_min, ra_max,
1942 dec_min, dec_max, field_margin);
1943 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1944 cpl_test_null(pErr);
1945
1946 pErr =
1947 hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, 0, ra_min, ra_max,
1948 dec_min, dec_max, field_margin);
1949 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1950 cpl_test_null(pErr);
1951
1952 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, 0, -1,
1953 ra_max, dec_min, dec_max, field_margin);
1954 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1955 cpl_test_null(pErr);
1956
1957 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, 0, ra_min,
1958 -1, dec_min, dec_max, field_margin);
1959 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1960 cpl_test_null(pErr);
1961
1962 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
1963 2, 1, dec_min, dec_max, field_margin);
1964 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1965 cpl_test_null(pErr);
1966
1967 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
1968 ra_min, ra_max, 2, 1, field_margin);
1969 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1970 cpl_test_null(pErr);
1971
1972 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
1973 ra_min, ra_max, dec_min, dec_max, -1);
1974 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1975 cpl_test_null(pErr);
1976
1977 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
1978 delta_lambda, -1, ra_max,
1979 dec_min, dec_max, lambda_min,
1980 lambda_max,field_margin);
1981 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1982 cpl_test_null(pErr);
1983
1984 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
1985 delta_lambda, ra_min, -1,
1986 dec_min, dec_max, lambda_min,
1987 lambda_max,field_margin);
1988 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1989 cpl_test_null(pErr);
1990
1991 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
1992 delta_lambda, 2, 1, dec_min,
1993 dec_max, lambda_min,
1994 lambda_max,field_margin);
1995 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
1996 cpl_test_null(pErr);
1997
1998 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
1999 delta_lambda, ra_min, ra_max,
2000 2, 1, lambda_min,
2001 lambda_max,field_margin);
2002 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2003 cpl_test_null(pErr);
2004
2005 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
2006 delta_lambda, ra_min, ra_max,
2007 dec_min, dec_max, -1,
2008 lambda_max,field_margin);
2009 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2010 cpl_test_null(pErr);
2011
2012 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
2013 delta_lambda, ra_min, ra_max,
2014 dec_min, dec_max, lambda_min,
2015 -1,field_margin);
2016 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2017 cpl_test_null(pErr);
2018
2019 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
2020 delta_lambda, ra_min, ra_max,
2021 dec_min, dec_max, 2,
2022 1,field_margin);
2023 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2024 cpl_test_null(pErr);
2025
2026 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
2027 delta_lambda, ra_min, ra_max,
2028 dec_min, dec_max, lambda_min,
2029 lambda_max,-1);
2030 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2031 cpl_test_null(pErr);
2032
2033 return cpl_error_get_code();
2034}
2035/*----------------------------------------------------------------------------*/
2040/*----------------------------------------------------------------------------*/
2041static cpl_error_code
2042test_resample_outgrid_param(void) {
2043 /* test invalid input */
2044 const double delta_ra = 0.1;
2045 const double delta_dec = 0.2;
2046 const double delta_lambda = 0.001;
2047
2048 hdrl_resample_outgrid_parameter * aParams_outputgrid =
2049 cpl_calloc (1, sizeof(hdrl_resample_outgrid_parameter));
2050
2051 aParams_outputgrid->delta_ra = delta_ra;
2052 aParams_outputgrid->delta_dec = delta_dec;
2053 aParams_outputgrid->delta_lambda = delta_lambda;
2054 aParams_outputgrid->wcs = hdrl_resample_util_crea_wcs_2d();
2055
2056 double ra_min = HDRL_RA_MIN; /* Minimal Right ascension [deg] */
2057 double ra_max = HDRL_RA_MIN; /* Maximal Right ascension [deg] */
2058 double dec_min = HDRL_DEC_MIN; /* Minimal Declination [deg] */
2059 double dec_max = HDRL_DEC_MIN; /* Maximal Declination [deg] */
2060 double lambda_min = HDRL_LAMBDA_MIN; /* Minimal wavelength [m] */
2061 double lambda_max = HDRL_LAMBDA_MIN; /* Maximal wavelength [m] */
2062 double field_margin = 5;
2063 aParams_outputgrid->dec_min = dec_min;
2064 aParams_outputgrid->dec_max = dec_max;
2065 aParams_outputgrid->ra_min = ra_min;
2066 aParams_outputgrid->ra_max = ra_max;
2067 aParams_outputgrid->lambda_min = lambda_min;
2068 aParams_outputgrid->lambda_max = lambda_max;
2069 hdrl_parameter* pErr = NULL;
2070 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
2071
2072 /* test invalid input */
2073 test_resample_invalid_input_outgrid_param(delta_ra, delta_dec, delta_lambda,
2074 ra_min, ra_max, dec_min, dec_max, lambda_min, lambda_max, field_margin);
2075
2076 /* test valid input */
2077 pErr = hdrl_resample_parameter_create_outgrid2D(delta_ra, delta_dec);
2078 cpl_test_error(CPL_ERROR_NONE);
2079 cpl_test_nonnull(pErr);
2080 cpl_test(hdrl_resample_parameter_outgrid_check(pErr));
2081 cpl_test_error(CPL_ERROR_NONE);
2082
2083 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_outgrid_verify(pErr));
2084 cpl_test_error(CPL_ERROR_NONE);
2085
2086 hdrl_parameter_delete(pErr);
2087 pErr = hdrl_resample_parameter_create_outgrid3D(delta_ra, delta_dec, delta_lambda);
2088 cpl_test_error(CPL_ERROR_NONE);
2089 cpl_test_nonnull(pErr);
2090 cpl_test(hdrl_resample_parameter_outgrid_check(pErr));
2091 cpl_test_error(CPL_ERROR_NONE);
2092
2093 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_outgrid_verify(pErr));
2094 cpl_test_error(CPL_ERROR_NONE);
2095
2096 hdrl_parameter_delete(pErr);
2097 pErr =
2098 hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
2099 ra_min, ra_max,
2100 dec_min, dec_max, field_margin);
2101 cpl_test_error(CPL_ERROR_NONE);
2102 cpl_test_nonnull(pErr);
2103 cpl_test(hdrl_resample_parameter_outgrid_check(pErr));
2104 cpl_test_error(CPL_ERROR_NONE);
2105
2106 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_outgrid_verify(pErr));
2107 cpl_test_error(CPL_ERROR_NONE);
2108
2109
2110 hdrl_parameter_delete(pErr);
2111 pErr = hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
2112 delta_lambda, ra_min, ra_max,
2113 dec_min, dec_max, lambda_min,
2114 lambda_max,field_margin);
2115 cpl_test_error(CPL_ERROR_NONE);
2116 cpl_test_nonnull(pErr);
2117
2118 cpl_test(hdrl_resample_parameter_outgrid_check(pErr));
2119 cpl_test_error(CPL_ERROR_NONE);
2120
2121 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_outgrid_verify(pErr));
2122 cpl_test_error(CPL_ERROR_NONE);
2123 hdrl_parameter_delete(pErr);
2124
2125 cpl_free(aParams_outputgrid);
2126
2127 return cpl_error_get_code();
2128}
2129/*----------------------------------------------------------------------------*/
2134/*----------------------------------------------------------------------------*/
2135static cpl_error_code
2136test_hdrl_resample_parameter_create_nearest(void) {
2137 /* test invalid input */
2138 hdrl_parameter *pErr;
2139
2140 pErr = hdrl_resample_parameter_create_nearest();
2141 cpl_test_error(CPL_ERROR_NONE);
2142 cpl_test_nonnull(pErr);
2143 hdrl_parameter_delete(pErr);
2144
2145 hdrl_parameter* p = hdrl_resample_parameter_create_nearest();
2146 cpl_test_nonnull(p);
2147
2148 cpl_test(hdrl_resample_parameter_method_check(p));
2149 cpl_test_error(CPL_ERROR_NONE);
2150
2151 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_method_verify(p));
2152 cpl_test_error(CPL_ERROR_NONE);
2153
2154 hdrl_parameter_delete(p);
2155
2156 return cpl_error_get_code();
2157}
2158/*----------------------------------------------------------------------------*/
2163/*----------------------------------------------------------------------------*/
2164static cpl_error_code
2165test_hdrl_resample_parameter_create_lanczos(void) {
2166 /* test invalid input */
2167 const double loop_distance = 2;
2168 const int kernel_size = 2;
2169 cpl_boolean use_errorweights= CPL_TRUE;
2170 hdrl_parameter *pErr;
2171 /* test invalid input */
2172 pErr = hdrl_resample_parameter_create_lanczos(1, use_errorweights, 0);
2173 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2174 pErr = hdrl_resample_parameter_create_lanczos(-1, use_errorweights,
2175 kernel_size);
2176 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2177 cpl_test_null(pErr);
2178 cpl_test_error(CPL_ERROR_NONE);
2179
2180 /* test valid input */
2181 hdrl_parameter* p = hdrl_resample_parameter_create_lanczos(loop_distance,
2182 use_errorweights,
2183 kernel_size);
2184 cpl_test_error(CPL_ERROR_NONE);
2185 cpl_test_nonnull(p);
2186
2187 cpl_test(hdrl_resample_parameter_method_check(p));
2188 cpl_test_error(CPL_ERROR_NONE);
2189
2190 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_method_verify(p));
2191 cpl_test_error(CPL_ERROR_NONE);
2192
2193 hdrl_parameter_delete(p);
2194
2195 return cpl_error_get_code();
2196}
2197/*----------------------------------------------------------------------------*/
2202/*----------------------------------------------------------------------------*/
2203static cpl_error_code
2204test_hdrl_resample_parameter_create_linear(void) {
2205 /* test invalid input */
2206 const double loop_distance = 2;
2207 cpl_boolean use_errorweights= CPL_TRUE;
2208
2209 hdrl_parameter *pErr;
2210 /* test invalid input */
2211 pErr = hdrl_resample_parameter_create_linear(-1, use_errorweights);
2212 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2213 cpl_test_null(pErr);
2214
2215 /* test valid input */
2216 hdrl_parameter* p = hdrl_resample_parameter_create_linear(loop_distance,
2217 use_errorweights);
2218 cpl_test_error(CPL_ERROR_NONE);
2219 cpl_test_nonnull(p);
2220
2221 cpl_test(hdrl_resample_parameter_method_check(p));
2222 cpl_test_error(CPL_ERROR_NONE);
2223
2224 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_method_verify(p));
2225 cpl_test_error(CPL_ERROR_NONE);
2226
2227 hdrl_parameter_delete(p);
2228
2229 return cpl_error_get_code();
2230}
2231/*----------------------------------------------------------------------------*/
2236/*----------------------------------------------------------------------------*/
2237static cpl_error_code
2238test_hdrl_resample_parameter_create_quadratic(void) {
2239 /* test invalid input */
2240 const double loop_distance = 2;
2241 cpl_boolean use_errorweights= CPL_TRUE;
2242 hdrl_parameter *pErr;
2243
2244 /* test invalid input */
2245 pErr = hdrl_resample_parameter_create_quadratic(-1, use_errorweights);
2246 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2247 cpl_test_null(pErr);
2248
2249 /* test valid input */
2250 hdrl_parameter* p = hdrl_resample_parameter_create_quadratic(loop_distance,
2251 use_errorweights);
2252 cpl_test_error(CPL_ERROR_NONE);
2253 cpl_test_nonnull(p);
2254
2255 cpl_test(hdrl_resample_parameter_method_check(p));
2256 cpl_test_error(CPL_ERROR_NONE);
2257
2258 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_method_verify(p));
2259 cpl_test_error(CPL_ERROR_NONE);
2260
2261 hdrl_parameter_delete(p);
2262
2263 return cpl_error_get_code();
2264}
2265/*----------------------------------------------------------------------------*/
2270/*----------------------------------------------------------------------------*/
2271static cpl_error_code
2272test_hdrl_resample_parameter_create_renka(void) {
2273 /* test invalid input */
2274 const double loop_distance = 2;
2275 const double critical_radius_renka = 3;
2276 cpl_boolean use_errorweights= CPL_TRUE;
2277 hdrl_parameter *pErr;
2278
2279 /* test invalid input */
2280 pErr = hdrl_resample_parameter_create_renka(-1, use_errorweights,
2281 critical_radius_renka);
2282 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2283 cpl_test_null(pErr);
2284
2285 pErr = hdrl_resample_parameter_create_renka(loop_distance,
2286 use_errorweights, -1);
2287 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2288 cpl_test_null(pErr);
2289
2290 /* test valid input */
2291 hdrl_parameter* p =
2292 hdrl_resample_parameter_create_renka(loop_distance, use_errorweights,
2293 critical_radius_renka);
2294
2295 cpl_test(hdrl_resample_parameter_method_check(p));
2296 cpl_test_error(CPL_ERROR_NONE);
2297
2298 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_method_verify(p));
2299 cpl_test_error(CPL_ERROR_NONE);
2300
2301 hdrl_parameter_delete(p);
2302
2303 return cpl_error_get_code();
2304}
2305static cpl_error_code
2306test_invalid_input_hdrl_resample_parameter_create_drizzle(
2307 const double loop_distance,
2308 const double pix_frac_drizzle_x,
2309 const double pix_frac_drizzle_y,
2310 const double pix_frac_drizzle_lambda,
2311 cpl_boolean use_errorweights
2312)
2313{
2314 hdrl_parameter *pErr;
2315
2316 /* test invalid input */
2317 pErr = hdrl_resample_parameter_create_drizzle(-1, use_errorweights,
2318 pix_frac_drizzle_x,
2319 pix_frac_drizzle_y,
2320 pix_frac_drizzle_lambda);
2321 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2322 cpl_test_null(pErr);
2323
2324 pErr = hdrl_resample_parameter_create_drizzle(loop_distance, use_errorweights,
2325 -1,
2326 pix_frac_drizzle_y,
2327 pix_frac_drizzle_lambda);
2328 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2329 cpl_test_null(pErr);
2330
2331 pErr = hdrl_resample_parameter_create_drizzle(loop_distance, use_errorweights,
2332 pix_frac_drizzle_x,
2333 -1,
2334 pix_frac_drizzle_lambda);
2335 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2336 cpl_test_null(pErr);
2337
2338 pErr = hdrl_resample_parameter_create_drizzle(loop_distance, use_errorweights,
2339 pix_frac_drizzle_x,
2340 pix_frac_drizzle_y,
2341 -1);
2342 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2343 cpl_test_null(pErr);
2344
2345 return cpl_error_get_code();
2346}
2347/*----------------------------------------------------------------------------*/
2352/*----------------------------------------------------------------------------*/
2353static cpl_error_code
2354test_hdrl_resample_parameter_create_drizzle(void) {
2355 /* test invalid input */
2356 const double loop_distance = 2;
2357 const double pix_frac_drizzle_x = 0.8;
2358 const double pix_frac_drizzle_y = 0.8;
2359 const double pix_frac_drizzle_lambda = 1;
2360 cpl_boolean use_errorweights= CPL_TRUE;
2361
2362 /* test valid input */
2363 test_invalid_input_hdrl_resample_parameter_create_drizzle(loop_distance,
2364 pix_frac_drizzle_x, pix_frac_drizzle_y, pix_frac_drizzle_lambda,
2365 use_errorweights);
2366
2367 /* test valid input */
2368 hdrl_parameter* p =
2369 hdrl_resample_parameter_create_drizzle(loop_distance, use_errorweights,
2370 pix_frac_drizzle_x,
2371 pix_frac_drizzle_y,
2372 pix_frac_drizzle_lambda);
2373
2374 cpl_test(hdrl_resample_parameter_method_check(p));
2375 cpl_test_error(CPL_ERROR_NONE);
2376
2377 cpl_test_eq(CPL_ERROR_NONE, hdrl_resample_parameter_method_verify(p));
2378 cpl_test_error(CPL_ERROR_NONE);
2379
2380 hdrl_parameter_delete(p);
2381
2382 return cpl_error_get_code();
2383}
2384/*----------------------------------------------------------------------------*/
2389/*----------------------------------------------------------------------------*/
2390static cpl_error_code
2391test_resample_parameters(void) {
2392
2393 test_resample_outgrid_param();
2394 cpl_test_error(CPL_ERROR_NONE);
2395
2396 test_hdrl_resample_parameter_create_nearest();
2397 cpl_test_error(CPL_ERROR_NONE);
2398
2399 test_hdrl_resample_parameter_create_lanczos();
2400 cpl_test_error(CPL_ERROR_NONE);
2401
2402 test_hdrl_resample_parameter_create_linear();
2403 cpl_test_error(CPL_ERROR_NONE);
2404
2405 test_hdrl_resample_parameter_create_quadratic();
2406 cpl_test_error(CPL_ERROR_NONE);
2407
2408 test_hdrl_resample_parameter_create_renka();
2409 cpl_test_error(CPL_ERROR_NONE);
2410
2411 test_hdrl_resample_parameter_create_drizzle();
2412 cpl_test_error(CPL_ERROR_NONE);
2413
2414 return cpl_error_get_code();
2415}
2416/*----------------------------------------------------------------------------*/
2421/*----------------------------------------------------------------------------*/
2422static cpl_error_code
2423test_hdrl_resample_parameter_create_outgrid2D(void){
2424 const double delta_ra = 0.1;
2425 const double delta_dec = 0.1;
2426 hdrl_parameter *pErr;
2427 /* test invalid input: nothing all values are allowed */
2428 /* test valid input */
2429 pErr = hdrl_resample_parameter_create_outgrid2D(delta_ra, delta_dec);
2430 cpl_test_error(CPL_ERROR_NONE);
2431 cpl_test_nonnull(pErr);
2432
2433 hdrl_resample_parameter_outgrid_verify(pErr);
2434 cpl_test_error(CPL_ERROR_NONE);
2435
2436 hdrl_resample_parameter_outgrid_check(pErr);
2437 cpl_test_error(CPL_ERROR_NONE);
2438 hdrl_parameter_delete(pErr);
2439 return cpl_error_get_code();
2440}
2441/*----------------------------------------------------------------------------*/
2446/*----------------------------------------------------------------------------*/
2447static cpl_error_code
2448test_hdrl_resample_parameter_create_outgrid3D(void){
2449 const double delta_ra = 0.1;
2450 const double delta_dec = 0.1;
2451 const double delta_lambda = 0.1;
2452 hdrl_parameter *pErr;
2453
2454 /* test invalid input: nothing all values are allowed */
2455
2456 /* test valid input */
2457 pErr = hdrl_resample_parameter_create_outgrid3D(delta_ra, delta_dec, delta_lambda);
2458 cpl_test_error(CPL_ERROR_NONE);
2459 cpl_test_nonnull(pErr);
2460
2461 hdrl_resample_parameter_outgrid_verify(pErr);
2462 cpl_test_error(CPL_ERROR_NONE);
2463
2464 hdrl_resample_parameter_outgrid_check(pErr);
2465 cpl_test_error(CPL_ERROR_NONE);
2466 hdrl_parameter_delete(pErr);
2467 return cpl_error_get_code();
2468}
2469
2470static cpl_error_code
2471test_invalid_input_hdrl_resample_parameter_create_outgrid2D_userdef(
2472 const double delta_ra,
2473 const double delta_dec,
2474 const double ra_min,
2475 const double ra_max,
2476 const double dec_min,
2477 const double dec_max,
2478 const double fieldmargin
2479){
2480
2481 hdrl_parameter *pErr;
2482
2483 /* test invalid input */
2484 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
2485 1, 0,
2486 dec_min, dec_max,
2487 fieldmargin);
2488 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2489 cpl_test_null(pErr);
2490
2491 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
2492 ra_min, ra_max,
2493 1, 0,
2494 fieldmargin);
2495 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2496 cpl_test_null(pErr);
2497
2498 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
2499 ra_min, ra_max,
2500 dec_min, dec_max,
2501 -1);
2502 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2503 cpl_test_null(pErr);
2504
2505 return cpl_error_get_code();
2506}
2507/*----------------------------------------------------------------------------*/
2512/*----------------------------------------------------------------------------*/
2513static cpl_error_code
2514test_hdrl_resample_parameter_create_outgrid2D_userdef(void){
2515 const double delta_ra = 0.1;
2516 const double delta_dec = 0.1;
2517 const double ra_min = 1;
2518 const double ra_max = 2;
2519 const double dec_min = 0;
2520 const double dec_max = 1;
2521 double fieldmargin = 5;
2522 hdrl_parameter *pErr;
2523
2524 /* test invalid input */
2525 test_invalid_input_hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra,
2526 delta_dec,
2527 ra_min, ra_max, dec_min, dec_max, fieldmargin);
2528 /* test valid input */
2529 pErr = hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
2530 ra_min, ra_max,
2531 dec_min, dec_max,
2532 fieldmargin);
2533 cpl_test_error(CPL_ERROR_NONE);
2534 cpl_test_nonnull(pErr);
2535
2536 hdrl_resample_parameter_outgrid_verify(pErr);
2537 cpl_test_error(CPL_ERROR_NONE);
2538
2539 hdrl_resample_parameter_outgrid_check(pErr);
2540 cpl_test_error(CPL_ERROR_NONE);
2541 hdrl_parameter_delete(pErr);
2542 return cpl_error_get_code();
2543}
2544static cpl_error_code
2545test_invalid_input_hdrl_resample_parameter_create_outgrid3D_userdef(
2546 const double delta_ra,
2547 const double delta_dec,
2548 const double delta_lambda,
2549 const double ra_min,
2550 const double ra_max,
2551 const double dec_min,
2552 const double dec_max,
2553 const double lambda_min,
2554 const double lambda_max,
2555 const double fieldmargin
2556){
2557 hdrl_parameter *pErr;
2558
2559 /* test invalid input */
2560 pErr =
2561 hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec, delta_lambda,
2562 2, 1,
2563 dec_min, dec_max,
2564 lambda_min, lambda_max,
2565 fieldmargin);
2566 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2567 cpl_test_null(pErr);
2568
2569 pErr =
2570 hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec, delta_lambda,
2571 ra_min, ra_max,
2572 2, 1,
2573 lambda_min, lambda_max,
2574 fieldmargin);
2575 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2576 cpl_test_null(pErr);
2577
2578 pErr =
2579 hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec, delta_lambda,
2580 ra_min, ra_max,
2581 dec_min, dec_max,
2582 550, 500,
2583 fieldmargin);
2584 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2585 cpl_test_null(pErr);
2586
2587 pErr =
2588 hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec, delta_lambda,
2589 ra_min, ra_max,
2590 dec_min, dec_max,
2591 lambda_min, lambda_max,
2592 -1);
2593 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2594 cpl_test_null(pErr);
2595
2596 return cpl_error_get_code();
2597}
2598/*----------------------------------------------------------------------------*/
2603/*----------------------------------------------------------------------------*/
2604static cpl_error_code
2605test_hdrl_resample_parameter_create_outgrid3D_userdef(void){
2606 const double delta_ra = 0.1;
2607 const double delta_dec = 0.1;
2608 const double delta_lambda = 0.1;
2609 const double ra_min = 1;
2610 const double ra_max = 2;
2611 const double dec_min = 0;
2612 const double dec_max = 1;
2613 const double lambda_min = 500;
2614 const double lambda_max = 550;
2615
2616 double fieldmargin = 5;
2617 hdrl_parameter *pErr;
2618
2619 /* test invalid input */
2620 test_invalid_input_hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra,
2621 delta_dec, delta_lambda, ra_min, ra_max, dec_min, dec_max, lambda_min,
2622 lambda_max, fieldmargin);
2623 /* test valid input */
2624 pErr =
2625 hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec, delta_lambda,
2626 ra_min, ra_max,
2627 dec_min, dec_max,
2628 lambda_min, lambda_max,
2629 fieldmargin);
2630 cpl_test_error(CPL_ERROR_NONE);
2631 cpl_test_nonnull(pErr);
2632
2633 hdrl_resample_parameter_outgrid_verify(pErr);
2634 cpl_test_error(CPL_ERROR_NONE);
2635
2636 hdrl_resample_parameter_outgrid_check(pErr);
2637 cpl_test_error(CPL_ERROR_NONE);
2638
2639 hdrl_parameter_delete(pErr);
2640 return cpl_error_get_code();
2641}
2642/*----------------------------------------------------------------------------*/
2647/*----------------------------------------------------------------------------*/
2648static cpl_error_code
2649test_hdrl_resample_wcs_print(void){
2650
2651 cpl_wcs* wcs = NULL;
2652
2653 /* verify invalid input */
2654 hdrl_resample_wcs_print(NULL);
2655 cpl_test_error(CPL_ERROR_NULL_INPUT);
2656
2657 /* verify valid input: 3D case */
2658 wcs = hdrl_resample_util_crea_wcs_3d();
2659 hdrl_resample_wcs_print(wcs);
2660 cpl_test_error(CPL_ERROR_NONE);
2661 cpl_wcs_delete(wcs);
2662
2663 /* verify valid input: 2D case */
2664 wcs = hdrl_resample_util_crea_wcs_2d();
2665 hdrl_resample_wcs_print(wcs);
2666 cpl_test_error(CPL_ERROR_NONE);
2667 cpl_wcs_delete(wcs);
2668
2669 return cpl_error_get_code();
2670}
2671/*----------------------------------------------------------------------------*/
2677/*----------------------------------------------------------------------------*/
2678static cpl_error_code
2679hdrl_resample_util_fill_outputgrid(hdrl_resample_outgrid_parameter * aParams_outputgrid)
2680{
2681 /*Assign the wcs*/
2682 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_2d();
2683 cpl_test_nonnull(wcs);
2684
2685 aParams_outputgrid->wcs = wcs;
2686
2687 /* Recalculate the limits if the user did not specify any */
2688 cpl_boolean recalc_limits = aParams_outputgrid->recalc_limits;
2689
2690 if (recalc_limits == CPL_TRUE) {
2691
2692 aParams_outputgrid->ra_min = HDRL_RA_MIN;
2693 aParams_outputgrid->ra_max = HDRL_RA_MAX;
2694 aParams_outputgrid->dec_min = HDRL_DEC_MIN;
2695 aParams_outputgrid->dec_max = HDRL_DEC_MAX;
2696 aParams_outputgrid->lambda_min = HDRL_LAMBDA_MIN;
2697 aParams_outputgrid->lambda_max = HDRL_LAMBDA_MAX;
2698 }
2699 cpl_test_error(CPL_ERROR_NONE);
2700 return cpl_error_get_code();
2701
2702}
2703/*----------------------------------------------------------------------------*/
2708/*----------------------------------------------------------------------------*/
2709static cpl_error_code
2710test_hdrl_resample_outgrid_parameter_print(void){
2711
2712 /* verify invalid input */
2713 hdrl_resample_outgrid_parameter_print(NULL);
2714 cpl_test_error(CPL_ERROR_NULL_INPUT);
2715
2716 hdrl_parameter * outputgrid =
2717 hdrl_resample_parameter_create_outgrid2D(1e-5, 1e-5);
2718
2719 hdrl_resample_outgrid_parameter *aParams_outputgrid =
2720 (hdrl_resample_outgrid_parameter *) outputgrid;
2721
2722 cpl_test_nonnull(aParams_outputgrid);
2723
2724 /*Assign the wcs*/
2725 hdrl_resample_util_fill_outputgrid(aParams_outputgrid);
2726 cpl_test_error(CPL_ERROR_NONE);
2727
2728 /* verify valid input */
2729 hdrl_resample_outgrid_parameter_print(aParams_outputgrid);
2730 cpl_test_error(CPL_ERROR_NONE);
2731
2732 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
2733 cpl_free(aParams_outputgrid);
2734
2735 return cpl_error_get_code();
2736}
2737/*----------------------------------------------------------------------------*/
2742/*----------------------------------------------------------------------------*/
2743static cpl_error_code
2744test_hdrl_resample_method_parameter_print(void){
2745
2746 /* verify invalid input */
2747 hdrl_resample_method_parameter_print(NULL);
2748 cpl_test_error(CPL_ERROR_NULL_INPUT);
2749
2750 hdrl_parameter * wrong = hdrl_resample_parameter_create_outgrid2D(0.1, 0.1);
2751 hdrl_resample_method_parameter * pwrong =
2752 (hdrl_resample_method_parameter *) wrong;
2753 cpl_test_nonnull(pwrong);
2754 hdrl_resample_method_parameter_print(pwrong);
2755 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
2756 cpl_free(pwrong);
2757
2758 /* verify valid input */
2759 hdrl_parameter * method = hdrl_resample_parameter_create_nearest();
2760 hdrl_resample_method_parameter * p =
2761 (hdrl_resample_method_parameter *) method;
2762 cpl_test_nonnull(p);
2763
2764 hdrl_resample_method_parameter_print(p);
2765 cpl_test_error(CPL_ERROR_NONE);
2766
2767 cpl_free(p);
2768 return cpl_error_get_code();
2769}
2770/*----------------------------------------------------------------------------*/
2775/*----------------------------------------------------------------------------*/
2776static cpl_error_code
2777test_hdrl_wcs_xy_to_radec(void) {
2778
2779 double x = 2;
2780 double y = 2;
2781 double ra = 2;
2782 double dec = 2;
2783
2784 /* verify invalid input */
2785 hdrl_wcs_xy_to_radec(NULL, x, y, &ra, &dec);
2786 cpl_test_error(CPL_ERROR_NULL_INPUT);
2787
2788 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_2d();
2789 cpl_test_nonnull(wcs);
2790
2791 hdrl_wcs_xy_to_radec(wcs, x, y, NULL, &dec);
2792 cpl_test_error(CPL_ERROR_NULL_INPUT);
2793
2794 hdrl_wcs_xy_to_radec(wcs, x, y, &ra, NULL);
2795 cpl_test_error(CPL_ERROR_NULL_INPUT);
2796
2797 /* verify valid input */
2798 hdrl_wcs_xy_to_radec(wcs, x, y, &ra, &dec);
2799 cpl_test_error(CPL_ERROR_NONE);
2800 cpl_msg_warning(cpl_func,"ra: %20.18g",ra);
2801 cpl_msg_warning(cpl_func,"dec: %20.18g",dec);
2802 cpl_test_abs(48.0716937492999961, ra, HDRL_EPS_TEST);
2803 cpl_test_abs(-20.6229937493000008, dec, HDRL_EPS_TEST);
2804
2805 cpl_wcs_delete(wcs);
2806 cpl_test_error(CPL_ERROR_NONE);
2807 return cpl_error_get_code();
2808}
2809/*----------------------------------------------------------------------------*/
2814/*----------------------------------------------------------------------------*/
2815static cpl_error_code
2816test_hdrl_resample_pfits_get(void) {
2817
2818 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_3d();
2819 cpl_test_nonnull(wcs);
2820 cpl_propertylist* header = cpl_propertylist_new();
2821 hdrl_wcs_to_propertylist(wcs, header, CPL_FALSE);
2822 cpl_test_nonnull(header);
2823 double crpix = 0;
2824 double crval = 0;
2825 double cd = 0;
2826
2827 /* test invalid input */
2828 crpix = hdrl_resample_pfits_get_crpix(NULL, 1);
2829 cpl_test_error(CPL_ERROR_NULL_INPUT);
2830 cpl_test_abs(crpix, 0, HDRL_EPS_TEST);
2831
2832 crval = hdrl_resample_pfits_get_crval(NULL, 1);
2833 cpl_test_error(CPL_ERROR_NULL_INPUT);
2834 cpl_test_abs(crval, 0, HDRL_EPS_TEST);
2835
2836 cd = hdrl_resample_pfits_get_cd(NULL, 1, 1);
2837 cpl_test_error(CPL_ERROR_NULL_INPUT);
2838 cpl_test_abs(cd, 0, HDRL_EPS_TEST);
2839
2840 /* test valid input */
2841 cpl_test_abs(HDRL_CRPIX1,hdrl_resample_pfits_get_crpix(header, 1), HDRL_EPS_TEST);
2842 cpl_test_abs(HDRL_CRPIX2,hdrl_resample_pfits_get_crpix(header, 2), HDRL_EPS_TEST);
2843 cpl_test_abs(HDRL_CRPIX3,hdrl_resample_pfits_get_crpix(header, 3), HDRL_EPS_TEST);
2844
2845 cpl_test_abs(HDRL_CRVAL1,hdrl_resample_pfits_get_crval(header, 1), HDRL_EPS_TEST);
2846 cpl_test_abs(HDRL_CRVAL2,hdrl_resample_pfits_get_crval(header, 2), HDRL_EPS_TEST);
2847 cpl_test_abs(HDRL_CRVAL3,hdrl_resample_pfits_get_crval(header, 3), HDRL_EPS_TEST);
2848
2849 cpl_test_abs(HDRL_CD11,hdrl_resample_pfits_get_cd(header, 1, 1), HDRL_EPS_TEST);
2850 cpl_test_abs(HDRL_CD12,hdrl_resample_pfits_get_cd(header, 1, 2), HDRL_EPS_TEST);
2851 cpl_test_abs(HDRL_CD21,hdrl_resample_pfits_get_cd(header, 2, 1), HDRL_EPS_TEST);
2852 cpl_test_abs(HDRL_CD22,hdrl_resample_pfits_get_cd(header, 2, 2), HDRL_EPS_TEST);
2853
2854 cpl_test_abs(HDRL_CD13,hdrl_resample_pfits_get_cd(header, 1, 3), HDRL_EPS_TEST);
2855 cpl_test_abs(HDRL_CD31,hdrl_resample_pfits_get_cd(header, 3, 1), HDRL_EPS_TEST);
2856 cpl_test_abs(HDRL_CD23,hdrl_resample_pfits_get_cd(header, 2, 3), HDRL_EPS_TEST);
2857 cpl_test_abs(HDRL_CD32,hdrl_resample_pfits_get_cd(header, 3, 1), HDRL_EPS_TEST);
2858 cpl_test_abs(HDRL_CD33,hdrl_resample_pfits_get_cd(header, 3, 3), HDRL_EPS_TEST);
2859
2860 cpl_propertylist_delete(header);
2861 cpl_wcs_delete(wcs);
2862 cpl_test_error(CPL_ERROR_NONE);
2863 return cpl_error_get_code();
2864}
2865/*----------------------------------------------------------------------------*/
2870/*----------------------------------------------------------------------------*/
2871static cpl_error_code
2872test_hdrl_resample_smallwcs_new(void){
2873
2874 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_3d();
2875 cpl_test_nonnull(wcs);
2876 cpl_propertylist* header = cpl_propertylist_new();
2877 hdrl_wcs_to_propertylist(wcs, header, CPL_FALSE);
2878 cpl_test_nonnull(header);
2879
2880 /* test invalid input */
2881 hdrl_resample_smallwcs* swcs = NULL;
2882 swcs = hdrl_resample_smallwcs_new(NULL);
2883 cpl_test_error(CPL_ERROR_NULL_INPUT);
2884 cpl_test_null(swcs);
2885
2886 /* test valid input */
2887 swcs = hdrl_resample_smallwcs_new(header);
2888 cpl_test_nonnull(swcs);
2889 cpl_test_error(CPL_ERROR_NONE);
2890 cpl_test_nonnull(swcs);
2891
2892 cpl_test_rel(swcs->cd11, HDRL_CD11, HDRL_EPS_TEST);
2893 cpl_test_rel(swcs->cd22, HDRL_CD22, HDRL_EPS_TEST);
2894 cpl_test_rel(swcs->cd12, HDRL_CD12, HDRL_EPS_TEST);
2895 cpl_test_rel(swcs->cd21, HDRL_CD21, HDRL_EPS_TEST);
2896 cpl_test_rel(swcs->crpix1, HDRL_CRPIX1, HDRL_EPS_TEST);
2897 cpl_test_rel(swcs->crpix2, HDRL_CRPIX2, HDRL_EPS_TEST);
2898 cpl_test_rel(swcs->crval1, HDRL_CRVAL1, HDRL_EPS_TEST);
2899 cpl_test_rel(swcs->crval2, HDRL_CRVAL2, HDRL_EPS_TEST);
2900
2901 cpl_free(swcs);
2902 cpl_propertylist_erase(header, "CRPIX1");
2903 swcs = hdrl_resample_smallwcs_new(header);
2904 cpl_test_nonnull(swcs);
2905 cpl_test_error(CPL_ERROR_NONE);
2906 cpl_test_nonnull(swcs);
2907 cpl_test_rel(swcs->crpix1, 0., HDRL_EPS_TEST);
2908
2909 cpl_free(swcs);
2910 cpl_propertylist_erase(header, "CD1_1");
2911 swcs = hdrl_resample_smallwcs_new(header);
2912 cpl_test_nonnull(swcs);
2913 cpl_test_error(CPL_ERROR_SINGULAR_MATRIX);
2914
2915 cpl_free(swcs);
2916 cpl_propertylist_erase_regexp(header, "CD?_?", 0);
2917 swcs = hdrl_resample_smallwcs_new(header);
2918 cpl_test_nonnull(swcs);
2919 cpl_test_error(CPL_ERROR_NONE);
2920 cpl_test_nonnull(swcs);
2921 cpl_test_rel(swcs->crpix1, 0., HDRL_EPS_TEST);
2922
2923 cpl_free(swcs);
2924 cpl_wcs_delete(wcs);
2925 cpl_propertylist_delete(header);
2926 cpl_test_error(CPL_ERROR_NONE);
2927 return cpl_error_get_code();
2928}
2929/*----------------------------------------------------------------------------*/
2934/*----------------------------------------------------------------------------*/
2935static hdrl_resample_pixgrid*
2936hdrl_resample_util_pixgrid_new(void) {
2937
2938 cpl_size aSizeX = 10;
2939 cpl_size aSizeY = 10;
2940 cpl_size aSizeZ = 10;
2941 unsigned short aNMaps = 10;
2942
2943 /* test valid input */
2944 hdrl_resample_pixgrid *aGrid =
2945 hdrl_resample_pixgrid_new(aSizeX, aSizeY, aSizeZ, aNMaps);
2946
2947 return aGrid;
2948}
2949/*----------------------------------------------------------------------------*/
2954/*----------------------------------------------------------------------------*/
2955static cpl_error_code
2956test_hdrl_resample_pixgrid_delete(void) {
2957
2958 /* test invalid input */
2959 hdrl_resample_pixgrid_delete(NULL);
2960 /* Nothing to check on error code as in this case function returns void */
2961
2962 /* test valid input */
2963 hdrl_resample_pixgrid *aGrid = hdrl_resample_util_pixgrid_new();
2964 cpl_test_nonnull(aGrid);
2965
2966 hdrl_resample_pixgrid_delete(aGrid);
2967 cpl_test_error(CPL_ERROR_NONE);
2968
2969 return cpl_error_get_code();
2970}
2971/*----------------------------------------------------------------------------*/
2976/*----------------------------------------------------------------------------*/
2977static cpl_table*
2978hdrl_resample_util_crea_pixel_table(void)
2979{
2980 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_3d();
2981 cpl_imagelist* ilist = NULL;
2982 cpl_imagelist* elist = NULL;
2983
2984 ilist = cpl_imagelist_new();
2985 elist = cpl_imagelist_new();
2986 cpl_image* simul = hdrl_resample_util_crea_image_dice_5(HDRL_SIZE_X, HDRL_SIZE_X,
2987 HDRL_FLUX_ADU);
2988 cpl_image* errs = cpl_image_duplicate(simul);
2989 cpl_image_power(errs, 0.5);
2990 cpl_imagelist_set(ilist, simul, 0);
2991 cpl_imagelist_set(elist, errs, 0);
2992 hdrl_imagelist* hlist = hdrl_resample_util_hdrl_imagelist_create(ilist,
2993 elist, NULL);
2994 cpl_table* pixel_table = hdrl_resample_imagelist_to_table(hlist, wcs);
2995
2996 cpl_imagelist_delete(ilist);
2997 cpl_imagelist_delete(elist);
2998 hdrl_imagelist_delete(hlist);
2999
3000 cpl_wcs_delete(wcs);
3001 cpl_test_error(CPL_ERROR_NONE);
3002
3003 return pixel_table;
3004}
3005/*----------------------------------------------------------------------------*/
3012/*----------------------------------------------------------------------------*/
3013
3014static cpl_error_code
3015test_hdrl_resample_compute_method(hdrl_resample_method aMethod,
3016 cpl_boolean use_errorweights) {
3017
3018 const double delta_ra = HDRL_CDELT1;
3019 const double delta_dec = HDRL_CDELT2;
3020 const double delta_lambda = HDRL_CDELT3;
3021 const double ramin = HDRL_RA_MIN;
3022 const double ramax = HDRL_RA_MAX;
3023 const double decmin = HDRL_DEC_MIN;
3024 const double decmax = HDRL_DEC_MAX;
3025 const double lambmin = HDRL_LAMBDA_MIN;
3026 const double lambmax = HDRL_LAMBDA_MAX;
3027
3028 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_3d();
3029
3030 cpl_table* pixel_table = hdrl_resample_util_crea_pixel_table();
3031
3032 hdrl_parameter * aParams_method =
3033 hdrl_resample_util_methodparam_create (
3034 LOOP_DISTANCE,
3035 RENKA_CRITICAL_RADIUS,
3036 LANCZOS_KERNEL_SIZE,
3037 DRIZZLE_DOWN_SCALING_FACTOR_X,
3038 DRIZZLE_DOWN_SCALING_FACTOR_Y,
3039 DRIZZLE_DOWN_SCALING_FACTOR_Z,
3040 aMethod, use_errorweights);
3041
3042 hdrl_parameter* aParams_outputgrid =
3043 hdrl_resample_parameter_create_outgrid3D_userdef(delta_ra, delta_dec,
3044 delta_lambda, ramin, ramax, decmin, decmax, lambmin, lambmax, 0.);
3045 /* we use 0 field margin to check later NAXISi values vs cube expected size */
3046
3047 /* Create resampling table starting from a data cube */
3048 int sx = HDRL_SIZE_X, sy = HDRL_SIZE_Z, sz = HDRL_SIZE_Z;
3049 cpl_image* ima = cpl_image_new(sx, sy, CPL_TYPE_INT);
3050 cpl_image_add_scalar(ima, 1);
3051 cpl_imagelist* iml = cpl_imagelist_new();
3052 for(int i = 0; i< sz;i++) {
3053 cpl_imagelist_set(iml, ima, i);
3054 }
3055 hdrl_imagelist* hlist = hdrl_resample_util_hdrl_imagelist_create(iml,
3056 NULL, NULL);
3057
3058 cpl_table * ResTable = hdrl_resample_imagelist_to_table(hlist, wcs);
3059 cpl_imagelist_delete(iml);
3060 hdrl_imagelist_delete(hlist);
3061
3062 cpl_test_error(CPL_ERROR_NONE);
3063 cpl_test_nonnull(ResTable);
3064 cpl_test_nonnull(aParams_method);
3065 cpl_test_nonnull(aParams_outputgrid);
3066 cpl_test_nonnull(wcs);
3067 hdrl_resample_result *cube =
3068 hdrl_resample_compute(ResTable, aParams_method, aParams_outputgrid,
3069 wcs);
3070 cpl_test_nonnull(cube);
3071 cpl_test_error(CPL_ERROR_NONE);
3072
3073 /* test valid input */
3074 /* test 3D case */
3075 hdrl_parameter_delete(aParams_outputgrid);
3076
3077 aParams_outputgrid =
3078 hdrl_resample_parameter_create_outgrid3D(delta_ra, delta_dec, delta_lambda);
3079
3080 /* test case: recalc_limits == CPL_TRUE */
3081 hdrl_resample_result_delete(cube);
3082 cube = hdrl_resample_compute(ResTable, aParams_method, aParams_outputgrid,
3083 wcs);
3084 cpl_test_nonnull(cube);
3085 cpl_test_error(CPL_ERROR_NONE);
3086
3087 /* test 2D case */
3088 hdrl_parameter_delete(aParams_outputgrid);
3089 aParams_outputgrid =
3090 hdrl_resample_parameter_create_outgrid2D_userdef(delta_ra, delta_dec,
3091 ramin, ramax,
3092 decmin, decmax,
3093 5.);
3094
3095 hdrl_resample_result_delete(cube);
3096 cube = hdrl_resample_compute(ResTable, aParams_method, aParams_outputgrid,
3097 wcs);
3098 cpl_test_nonnull(cube);
3099 cpl_test_error(CPL_ERROR_NONE);
3100
3101 hdrl_parameter_delete(aParams_outputgrid);
3102 hdrl_parameter_delete(aParams_method);
3103 cpl_table_delete(ResTable);
3104 cpl_wcs_delete(wcs);
3105 cpl_table_delete(pixel_table);
3106 hdrl_resample_result_delete(cube);
3107
3108 return cpl_error_get_code();
3109}
3110/*----------------------------------------------------------------------------*/
3115/*----------------------------------------------------------------------------*/
3116static cpl_error_code
3117test_hdrl_resample_compute(void) {
3118
3119 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_NEAREST, CPL_TRUE);
3120 cpl_test_error(CPL_ERROR_NONE);
3121 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_RENKA, CPL_TRUE);
3122 cpl_test_error(CPL_ERROR_NONE);
3123 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_LINEAR, CPL_TRUE);
3124 cpl_test_error(CPL_ERROR_NONE);
3125 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_QUADRATIC, CPL_TRUE);
3126 cpl_test_error(CPL_ERROR_NONE);
3127 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_DRIZZLE, CPL_TRUE);
3128 cpl_test_error(CPL_ERROR_NONE);
3129 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_LANCZOS, CPL_TRUE);
3130 cpl_test_error(CPL_ERROR_NONE);
3131 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_NONE, CPL_TRUE);
3132 cpl_test_error(CPL_ERROR_NONE);
3133 test_hdrl_resample_compute_method(HDRL_RESAMPLE_METHOD_LINEAR, CPL_FALSE);
3134 cpl_test_error(CPL_ERROR_NONE);
3135
3136 return cpl_error_get_code();
3137}
3138/*----------------------------------------------------------------------------*/
3143/*----------------------------------------------------------------------------*/
3144static cpl_error_code
3145test_hdrl_resample_wcs_projplane_from_celestial(void) {
3146
3147 hdrl_parameter * outputgrid =
3148 hdrl_resample_parameter_create_outgrid2D(1e-5, 1e-5);
3149
3150 hdrl_resample_outgrid_parameter *pogrid =
3151 (hdrl_resample_outgrid_parameter *) outputgrid;
3152
3153 double ra = HDRL_RA;
3154 double dec = HDRL_DEC;
3155 double x_out = HDRL_RA;
3156 double y_out = HDRL_DEC;
3157
3158 /* test invalid input */
3159 hdrl_resample_wcs_projplane_from_celestial(NULL, ra, dec, &x_out, &y_out);
3160 cpl_test_error(CPL_ERROR_NULL_INPUT);
3161
3162 /* test valid input */
3163 hdrl_resample_util_fill_outputgrid(pogrid);
3164 hdrl_resample_wcs_projplane_from_celestial(pogrid, ra, dec, &x_out, &y_out);
3165 cpl_test_error(CPL_ERROR_NONE);
3166 cpl_test_abs(x_out, -0.000561558354797, HDRL_DELTA_COMPARE_VALUE_ABS);
3167 cpl_test_abs(y_out, 0.000899999036676, HDRL_DELTA_COMPARE_VALUE_ABS);
3168
3169 cpl_wcs_delete((cpl_wcs *)pogrid->wcs);
3170 cpl_free(pogrid);
3171 cpl_test_error(CPL_ERROR_NONE);
3172 return cpl_error_get_code();
3173
3174}
3175/*----------------------------------------------------------------------------*/
3180/*----------------------------------------------------------------------------*/
3181static cpl_error_code
3182test_hdrl_resample_wcs_pixel_from_celestial_fast(void){
3183
3184 double ra = HDRL_RA;
3185 double dec = HDRL_DEC;
3186 double x_out = HDRL_RA;
3187 double y_out = HDRL_DEC;
3188 /* test invalid input */
3189 hdrl_resample_wcs_pixel_from_celestial_fast(NULL, ra, dec, &x_out, &y_out);
3190 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_3d();
3191 cpl_test_nonnull(wcs);
3192 cpl_propertylist* header = cpl_propertylist_new();
3193 hdrl_wcs_to_propertylist(wcs, header, CPL_FALSE);
3194 cpl_test_nonnull(header);
3195
3196 hdrl_resample_smallwcs * swcs = hdrl_resample_smallwcs_new(header);
3197 cpl_test_nonnull(swcs);
3198
3199 hdrl_resample_wcs_pixel_from_celestial_fast(swcs, ra, dec, &x_out, &y_out);
3200 cpl_test_error(CPL_ERROR_NONE);
3201
3202 /*
3203 This function uses many trigonometric functions and thus accumulates differences between
3204 e.g. linux 64 bit, linux 32 bit and macosx
3205 */
3206 cpl_test_rel(x_out, -163.82544787203, 1.e-11);
3207 cpl_test_rel(y_out, 1518.66596649235, 1.e-11);
3208
3209 cpl_free(swcs);
3210 cpl_wcs_delete(wcs);
3211 cpl_propertylist_delete(header);
3212 cpl_test_error(CPL_ERROR_NONE);
3213 return cpl_error_get_code();
3214}
3215
3216/*----------------------------------------------------------------------------*/
3221/*----------------------------------------------------------------------------*/
3222static cpl_error_code
3223test_hdrl_resample_compute_size(void){
3224
3225 hdrl_parameter * outputgrid =
3226 hdrl_resample_parameter_create_outgrid2D(1e-5, 1e-5);
3227
3228 hdrl_resample_outgrid_parameter *pogrid =
3229 (hdrl_resample_outgrid_parameter *) outputgrid;
3230
3231 int x, y, z;
3232
3233 /* test invalid input */
3234 hdrl_resample_compute_size(NULL, &x, &y, &z);
3235 cpl_test_error(CPL_ERROR_NULL_INPUT);
3236
3237 /* test valid input */
3238 hdrl_resample_util_fill_outputgrid(pogrid);
3239 hdrl_resample_compute_size(pogrid, &x, &y, &z);
3240 cpl_test_error(CPL_ERROR_NONE);
3241
3242 cpl_test_eq(x, 225); /* values determined on fc32 */
3243 cpl_test_eq(y, 229); /* values determined on fc32 */
3244 cpl_test_eq(z, 2); /* values determined on fc32 */
3245
3246 cpl_wcs_delete((cpl_wcs *)pogrid->wcs);
3247 cpl_free(pogrid);
3248 cpl_test_error(CPL_ERROR_NONE);
3249 return cpl_error_get_code();
3250}
3251/*----------------------------------------------------------------------------*/
3280/*----------------------------------------------------------------------------*/
3281static cpl_error_code
3282test_hdrl_resample_compute_size_rotated_wcs_bug(void)
3283{
3284 /* Problematic WCS CD matrix values from real dataset */
3285 const double cd11 = -3.80403139678984E-08;
3286 const double cd12 = -6.99483821074385E-05; /* Large off-diagonal */
3287 const double cd21 = -6.99981756546815E-05; /* Large off-diagonal */
3288 const double cd22 = 4.28668917003206E-08;
3289
3290 /* Incorrect delta values using only diagonal elements (the bug) */
3291 const double incorrect_delta_ra = fabs(cd11); /* 3.804e-08 deg */
3292 const double incorrect_delta_dec = fabs(cd22); /* 4.287e-08 deg */
3293
3294 /* Correct delta values accounting for rotation */
3295 const double correct_delta_ra = sqrt(cd11 * cd11 + cd12 * cd12); /* ~6.995e-05 deg */
3296 const double correct_delta_dec = sqrt(cd21 * cd21 + cd22 * cd22); /* ~6.999e-05 deg */
3297
3298 /* Create WCS with rotated CD matrix */
3299 cpl_propertylist* plist = cpl_propertylist_new();
3300 cpl_propertylist_append_int(plist, "NAXIS", 2);
3301 cpl_propertylist_append_int(plist, "NAXIS1", 100);
3302 cpl_propertylist_append_int(plist, "NAXIS2", 100);
3303 cpl_propertylist_append_double(plist, "CRPIX1", 50.0);
3304 cpl_propertylist_append_double(plist, "CRPIX2", 50.0);
3305 cpl_propertylist_append_double(plist, "CRVAL1", 48.0706);
3306 cpl_propertylist_append_double(plist, "CRVAL2", -20.6219);
3307 cpl_propertylist_append_string(plist, "CTYPE1", "RA---TAN");
3308 cpl_propertylist_append_string(plist, "CTYPE2", "DEC--TAN");
3309 cpl_propertylist_append_string(plist, "CUNIT1", "deg");
3310 cpl_propertylist_append_string(plist, "CUNIT2", "deg");
3311
3312 /* Set CD matrix with rotation */
3313 cpl_propertylist_append_double(plist, "CD1_1", cd11);
3314 cpl_propertylist_append_double(plist, "CD1_2", cd12);
3315 cpl_propertylist_append_double(plist, "CD2_1", cd21);
3316 cpl_propertylist_append_double(plist, "CD2_2", cd22);
3317
3318 /* Save and reload to ensure proper FITS header format */
3319 cpl_image* ima = cpl_image_new(100, 100, CPL_TYPE_INT);
3320 cpl_image_add_scalar(ima, 1);
3321 cpl_image_save(ima, "test_rotated_wcs.fits", CPL_TYPE_INT, plist, CPL_IO_DEFAULT);
3322 cpl_image_delete(ima);
3323 cpl_propertylist_delete(plist);
3324 plist = cpl_propertylist_load("test_rotated_wcs.fits", 0);
3325
3326 cpl_wcs* wcs = cpl_wcs_new_from_propertylist(plist);
3327 cpl_test_nonnull(wcs);
3328 cpl_test_error(CPL_ERROR_NONE);
3329 cpl_propertylist_delete(plist);
3330
3331 /* Create outputgrid parameter with INCORRECT delta values (simulating the bug) */
3332 hdrl_parameter* outputgrid =
3333 hdrl_resample_parameter_create_outgrid2D(incorrect_delta_ra, incorrect_delta_dec);
3334 cpl_test_nonnull(outputgrid);
3335 cpl_test_error(CPL_ERROR_NONE);
3336
3337 hdrl_resample_outgrid_parameter* pogrid =
3338 (hdrl_resample_outgrid_parameter*) outputgrid;
3339
3340 /* Set WCS and coordinate limits */
3341 pogrid->wcs = wcs;
3342 pogrid->ra_min = 48.069416667;
3343 pogrid->ra_max = 48.0718125;
3344 pogrid->dec_min = -20.6229925;
3345 pogrid->dec_max = -20.620708611;
3346
3347 int x, y, z;
3348
3349 /* Call hdrl_resample_compute_size with incorrect delta values */
3350 /* This should expose the bug by computing unreasonably large grid dimensions */
3351 hdrl_resample_compute_size(pogrid, &x, &y, &z);
3352 cpl_test_error(CPL_ERROR_NONE);
3353
3354 /* Calculate expected grid size with correct delta values */
3355 int expected_x = (int)(fabs(pogrid->ra_max - pogrid->ra_min) / correct_delta_ra) + 1;
3356 int expected_y = (int)(fabs(pogrid->dec_max - pogrid->dec_min) / correct_delta_dec) + 1;
3357
3358 /* The bug manifests as extremely large grid dimensions */
3359 /* With incorrect delta values (~3.8e-08), the grid size would be huge */
3360 /* A reasonable grid for this field should be around 1000-2000 pixels */
3361 /* With the bug, it would be orders of magnitude larger */
3362 const int MAX_REASONABLE_SIZE = 10000; /* Reasonable upper limit */
3363
3364 cpl_msg_info(cpl_func, "PIPE-12792 Bug Test Results:");
3365 cpl_msg_info(cpl_func, " Incorrect delta_ra: %g deg (diagonal only)", incorrect_delta_ra);
3366 cpl_msg_info(cpl_func, " Incorrect delta_dec: %g deg (diagonal only)", incorrect_delta_dec);
3367 cpl_msg_info(cpl_func, " Correct delta_ra: %g deg (with rotation)", correct_delta_ra);
3368 cpl_msg_info(cpl_func, " Correct delta_dec: %g deg (with rotation)", correct_delta_dec);
3369 cpl_msg_info(cpl_func, " Computed grid size: %d x %d x %d", x, y, z);
3370 cpl_msg_info(cpl_func, " Expected grid size with correct delta values: %d x %d",
3371 expected_x, expected_y);
3372
3373 /* This test FAILS when the bug exists (grid dimensions unreasonably large) */
3374 /* After the fix is applied, the test should PASS (reasonable grid dimensions) */
3375 /* The test verifies that grid dimensions are reasonable and match expected values */
3376
3377 /* Check that grid dimensions are not unreasonably large (bug check) */
3378 if (x > MAX_REASONABLE_SIZE || y > MAX_REASONABLE_SIZE) {
3379 cpl_msg_warning(cpl_func,
3380 "BUG DETECTED: Grid size (%d x %d) is unreasonably large! "
3381 "This indicates the pixel scale bug (PIPE-12792). "
3382 "Expected size with correct delta values: ~%d x %d",
3383 x, y, expected_x, expected_y);
3384 }
3385
3386 /* Test FAILS if grid dimensions are unreasonably large (bug exists) */
3387 /* Test PASSES if grid dimensions are reasonable (bug fixed) */
3388 cpl_test(x <= MAX_REASONABLE_SIZE && y <= MAX_REASONABLE_SIZE);
3389
3390 /* Additionally verify that grid dimensions match expected reasonable values */
3391 /* Allow some tolerance for rounding differences */
3392 const int TOLERANCE = 5; /* Allow ±5 pixels tolerance */
3393 cpl_test_abs(x, expected_x, TOLERANCE);
3394 cpl_test_abs(y, expected_y, TOLERANCE);
3395
3396 cpl_wcs_delete(wcs);
3397 cpl_free(pogrid);
3398 cpl_test_error(CPL_ERROR_NONE);
3399 return cpl_error_get_code();
3400}
3401/*----------------------------------------------------------------------------*/
3406/*----------------------------------------------------------------------------*/
3407static cpl_error_code
3408test_hdrl_resample_pixgrid_add(void){
3409
3410 hdrl_resample_pixgrid *aGrid = hdrl_resample_util_pixgrid_new();
3411 cpl_test_nonnull(aGrid);
3412
3413 cpl_size aIndex = 1;
3414 cpl_size aRow = 1;
3415 unsigned short aXIdx = 1;
3416
3417 /* test invalid input */
3418 hdrl_resample_pixgrid_add(NULL, aIndex, aRow, aXIdx);
3419 hdrl_resample_pixgrid_add(aGrid, -1, aRow, aXIdx);
3420
3421 /* test valid input */
3422 hdrl_resample_pixgrid_add(aGrid, aIndex, aRow, aXIdx);
3423 /* As hdrl_resample_pixgrid_add returns void we just test error code */
3424 cpl_test_error(CPL_ERROR_NONE);
3425
3426 hdrl_resample_pixgrid_delete(aGrid);
3427 cpl_test_error(CPL_ERROR_NONE);
3428
3429 return cpl_error_get_code();
3430}
3431/*----------------------------------------------------------------------------*/
3436/*----------------------------------------------------------------------------*/
3437static cpl_error_code
3438test_hdrl_resample_pixgrid_get_count(void){
3439 cpl_size aIndex = 1;
3440 /* test invalid input */
3441 hdrl_resample_pixgrid_get_count(NULL, aIndex);
3442
3443 hdrl_resample_pixgrid *aGrid = hdrl_resample_util_pixgrid_new();
3444 cpl_test_nonnull(aGrid);
3445
3446 hdrl_resample_pixgrid_get_count(aGrid, -1);
3447
3448 /* test valid input */
3449 cpl_size n = hdrl_resample_pixgrid_get_count(aGrid, aIndex);
3450
3451 cpl_test_rel(n, 0, HDRL_EPS_TEST);
3452
3453 hdrl_resample_pixgrid_delete(aGrid);
3454 cpl_test_error(CPL_ERROR_NONE);
3455 return cpl_error_get_code();
3456
3457}
3458/*----------------------------------------------------------------------------*/
3463/*----------------------------------------------------------------------------*/
3464static cpl_error_code
3465test_hdrl_resample_pixgrid_get_index(void)
3466{
3467 cpl_size x = 1;
3468 cpl_size y = 1;
3469 cpl_size z = 1;
3470 cpl_boolean aAllowOutside = CPL_TRUE;
3471
3472 hdrl_resample_pixgrid *aGrid = hdrl_resample_util_pixgrid_new();
3473 cpl_test_nonnull(aGrid);
3474 cpl_msg_warning(cpl_func,"x: %lld",aGrid->nx);
3475 cpl_msg_warning(cpl_func,"y: %lld",aGrid->ny);
3476 cpl_msg_warning(cpl_func,"z: %lld",aGrid->nz);
3477 cpl_size index = 0;
3478
3479 /* test invalid input */
3480 index = hdrl_resample_pixgrid_get_index(NULL, x, y, z, aAllowOutside);
3481 cpl_test_eq(index, -1);
3482 index = hdrl_resample_pixgrid_get_index(aGrid, -1, y, z, CPL_FALSE);
3483 cpl_test_eq(index, -1);
3484 index = hdrl_resample_pixgrid_get_index(aGrid, x, -1, z, CPL_FALSE);
3485 cpl_test_eq(index, -1);
3486 index = hdrl_resample_pixgrid_get_index(aGrid, x, y, -1, CPL_FALSE);
3487 cpl_test_eq(index, -1);
3488
3489 index = hdrl_resample_pixgrid_get_index(aGrid, aGrid->nx+1, y, z, CPL_FALSE);
3490 cpl_test_eq(index, -1);
3491 index = hdrl_resample_pixgrid_get_index(aGrid, x, aGrid->ny+1, z, CPL_FALSE);
3492 cpl_test_eq(index, -1);
3493 index = hdrl_resample_pixgrid_get_index(aGrid, x, y, aGrid->nz+1, CPL_FALSE);
3494 cpl_test_eq(index, -1);
3495
3496 index = hdrl_resample_pixgrid_get_index(aGrid, aGrid->nx + 1, y, z, CPL_TRUE);
3497 cpl_test_eq(index, 119);
3498 index = hdrl_resample_pixgrid_get_index(aGrid, x, aGrid->ny + 1, z, CPL_TRUE);
3499 cpl_test_eq(index, 191);
3500 index = hdrl_resample_pixgrid_get_index(aGrid, x, y, aGrid->nz + 1, CPL_TRUE);
3501 cpl_test_eq(index, 911);
3502
3503 /* test valid input */
3504 cpl_size index_res = 111;
3505 cpl_msg_warning(cpl_func,"index: %lld",
3506 hdrl_resample_pixgrid_get_index(aGrid, x, y, z, aAllowOutside));
3507 index = hdrl_resample_pixgrid_get_index(aGrid, x, y, z, aAllowOutside);
3508 cpl_test_eq(index_res, index);
3509
3510 hdrl_resample_pixgrid_delete(aGrid);
3511 cpl_test_error(CPL_ERROR_NONE);
3512 return cpl_error_get_code();
3513}
3514/*----------------------------------------------------------------------------*/
3519/*----------------------------------------------------------------------------*/
3520static cpl_error_code
3521test_hdrl_resample_pixgrid_new(void) {
3522 cpl_size x = 1;
3523 cpl_size y = 1;
3524 cpl_size z = 1;
3525 unsigned short aNMaps = 1;
3526 hdrl_resample_pixgrid * pGrid = NULL;
3527
3528 /* test invalid input */
3529 pGrid = hdrl_resample_pixgrid_new(-1, y, z, aNMaps);
3530 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3531 cpl_test_null(pGrid);
3532
3533 pGrid = hdrl_resample_pixgrid_new(x, -1, z, aNMaps);
3534 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3535 cpl_test_null(pGrid);
3536
3537 pGrid = hdrl_resample_pixgrid_new(x, y, -1, aNMaps);
3538 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3539 cpl_test_null(pGrid);
3540
3541 pGrid = hdrl_resample_pixgrid_new(x, y, z, 0);
3542 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3543 cpl_test_null(pGrid);
3544
3545 /* test valid input */
3546 pGrid = hdrl_resample_pixgrid_new(x, y, z, aNMaps);
3547 cpl_test_error(CPL_ERROR_NONE);
3548 cpl_test_nonnull(pGrid);
3549
3550 hdrl_resample_pixgrid_delete(pGrid);
3551 return cpl_error_get_code();
3552}
3553/*----------------------------------------------------------------------------*/
3558/*----------------------------------------------------------------------------*/
3559static cpl_error_code
3560test_hdrl_resample_pixgrid_create(void) {
3561
3562 cpl_size x = 1;
3563 cpl_size y = 1;
3564 cpl_size z = 1;
3565 cpl_table* pixel_table = NULL;
3566 cpl_wcs* wcs = hdrl_resample_util_crea_wcs_3d();
3567 cpl_test_nonnull(wcs);
3568 cpl_propertylist* header = cpl_propertylist_new();
3569 hdrl_wcs_to_propertylist(wcs, header, CPL_FALSE);
3570 hdrl_resample_pixgrid * pGrid= NULL;
3571
3572 /* test invalid input */
3573 pixel_table = cpl_table_new(0);
3574 pGrid = hdrl_resample_pixgrid_create(pixel_table, header, x, y, z);
3575 cpl_test_error(CPL_ERROR_NULL_INPUT);
3576 cpl_test_null(pGrid);
3577 cpl_table_delete(pixel_table);
3578
3579 pixel_table = cpl_table_new(1);
3580 pGrid = hdrl_resample_pixgrid_create(pixel_table, header, x, y, z);
3581 cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
3582 cpl_test_null(pGrid);
3583 cpl_table_delete(pixel_table);
3584
3585 pixel_table = hdrl_resample_util_crea_pixel_table();
3586
3587 /* test invalid input */
3588 pGrid = hdrl_resample_pixgrid_create(pixel_table, NULL, x, y, z);
3589 cpl_test_error(CPL_ERROR_NULL_INPUT);
3590 cpl_test_null(pGrid);
3591
3592 pGrid = hdrl_resample_pixgrid_create(pixel_table, header, -1, y, z);
3593 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3594 cpl_test_null(pGrid);
3595
3596 pGrid = hdrl_resample_pixgrid_create(pixel_table, header, x, -1, z);
3597 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3598 cpl_test_null(pGrid);
3599
3600 pGrid = hdrl_resample_pixgrid_create(pixel_table, header, x, y, -1);
3601 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3602 cpl_test_null(pGrid);
3603
3604 /* test valid input */
3605 pGrid = hdrl_resample_pixgrid_create(pixel_table, header, x, y, z);
3606 cpl_test_error(CPL_ERROR_NONE);
3607 cpl_test_nonnull(pGrid);
3608 hdrl_resample_pixgrid_delete(pGrid);
3609
3610 cpl_table_unselect_row(pixel_table, 1);
3611 pGrid = hdrl_resample_pixgrid_create(pixel_table, header, x, y, z);
3612 cpl_test_error(CPL_ERROR_NONE);
3613 cpl_test_nonnull(pGrid);
3614
3615 cpl_wcs_delete(wcs);
3616 cpl_table_delete(pixel_table);
3617 cpl_propertylist_delete(header);
3618 hdrl_resample_pixgrid_delete(pGrid);
3619 return cpl_error_get_code();
3620}
3621/*----------------------------------------------------------------------------*/
3626/*----------------------------------------------------------------------------*/
3627static cpl_error_code
3628test_hdrl_resample_wcs_get_scales(void) {
3629
3630 double aXScale = 1;
3631 double aYScale = 1;
3632 hdrl_parameter * outputgrid =
3633 hdrl_resample_parameter_create_outgrid2D(1e-5, 1e-5);
3634
3635 hdrl_resample_outgrid_parameter *pogrid =
3636 (hdrl_resample_outgrid_parameter *) outputgrid;
3637 hdrl_resample_util_fill_outputgrid(pogrid);
3638 hdrl_parameter * method = hdrl_resample_parameter_create_nearest();
3639
3640 hdrl_resample_method_parameter * p =
3641 (hdrl_resample_method_parameter *) method;
3642 cpl_test_nonnull(p);
3643
3644 /* test invalid input */
3645 hdrl_resample_wcs_get_scales(NULL, &aXScale, &aYScale);
3646 cpl_test_error(CPL_ERROR_NULL_INPUT);
3647
3648 hdrl_resample_wcs_get_scales(pogrid, NULL, &aYScale);
3649 cpl_test_error(CPL_ERROR_NULL_INPUT);
3650
3651 hdrl_resample_wcs_get_scales(pogrid, &aXScale, NULL);
3652 cpl_test_error(CPL_ERROR_NULL_INPUT);
3653
3654 /* test valid input */
3655 hdrl_resample_wcs_get_scales(pogrid, &aXScale, &aYScale);
3656 cpl_msg_warning(cpl_func,"x: %16.10g",aXScale);
3657 cpl_msg_warning(cpl_func,"y: %16.10g",aYScale);
3658 cpl_test_error(CPL_ERROR_NONE);
3659 /* clean memory */
3660 cpl_wcs_delete((cpl_wcs *)pogrid->wcs);
3661 hdrl_parameter_delete(outputgrid);
3662 hdrl_parameter_delete(method);
3663 return cpl_error_get_code();
3664}
3665/*----------------------------------------------------------------------------*/
3670/*----------------------------------------------------------------------------*/
3671static cpl_error_code
3672test_hdrl_resample_create_table(void){
3673
3674 const cpl_size size = 10;
3675 cpl_table* tab = NULL;
3676
3677 /* test invalid input */
3678 hdrl_resample_create_table(NULL, size);
3679 cpl_test_error(CPL_ERROR_NULL_INPUT);
3680
3681 hdrl_resample_create_table(&tab, -1);
3682 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3683
3684 /* test valid input */
3685 hdrl_resample_create_table(&tab, size);
3686 cpl_test_error(CPL_ERROR_NONE);
3687
3688 cpl_table_delete(tab);
3689 return cpl_error_get_code();
3690}
3691/*----------------------------------------------------------------------------*/
3696/*----------------------------------------------------------------------------*/
3697static cpl_error_code
3698test_hdrl_resample_pixgrid_get_rows(void) {
3699
3700 cpl_size aIndex = 1;
3701 hdrl_resample_pixgrid *aGrid = NULL;
3702 cpl_size x = 2;
3703 cpl_size y = 1;
3704 cpl_size z = 1;
3705 unsigned short aNMaps = 1;
3706
3707 /* test invalid input */
3708 hdrl_resample_pixgrid_get_rows(NULL, aIndex);
3709 cpl_test_error(CPL_ERROR_NULL_INPUT);
3710
3711 aGrid = hdrl_resample_pixgrid_new(x, y, z, aNMaps);
3712 cpl_test_error(CPL_ERROR_NONE);
3713 cpl_test_nonnull(aGrid);
3714
3715 hdrl_resample_pixgrid_get_rows(aGrid, -1);
3716 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
3717
3718 /* test valid input */
3719 hdrl_resample_pixgrid_get_rows(aGrid, aIndex);
3720 cpl_test_error(CPL_ERROR_NONE);
3721 hdrl_resample_pixgrid_delete(aGrid);
3722
3723 return cpl_error_get_code();
3724}
3725/*----------------------------------------------------------------------------*/
3730/*----------------------------------------------------------------------------*/
3731static cpl_error_code
3732test_hdrl_resample_cube_nearest(void){
3733
3734 cpl_table *ResTable = NULL;
3735 const cpl_size size = 10;
3736 cpl_size x = 2;
3737 cpl_size y = 2;
3738 cpl_size z = 2;
3739 unsigned short aNMaps = 1;
3740
3741 /* prepare re-sample table */
3742 hdrl_resample_create_table(&ResTable, size);
3743
3744 /* prepare a Grid */
3745 hdrl_resample_pixgrid *aGrid = hdrl_resample_pixgrid_new(x, y, z, aNMaps);
3746
3747 /* prepare output-grid parameter */
3748 hdrl_parameter * outputgrid =
3749 hdrl_resample_parameter_create_outgrid3D_userdef(HDRL_CDELT1,
3750 HDRL_CDELT2, HDRL_CDELT3,
3751 HDRL_RA_MIN, HDRL_RA_MAX, HDRL_DEC_MIN, HDRL_DEC_MAX,
3752 HDRL_LAMBDA_MIN, HDRL_LAMBDA_MAX,5.);
3753
3754 /* creates the cube */
3755 hdrl_resample_method aMethod = HDRL_RESAMPLE_METHOD_DRIZZLE;
3756 cpl_boolean use_errorweights = CPL_TRUE;
3757 hdrl_parameter * aParams_method =
3758 hdrl_resample_util_methodparam_create (
3759 LOOP_DISTANCE,
3760 RENKA_CRITICAL_RADIUS,
3761 LANCZOS_KERNEL_SIZE,
3762 DRIZZLE_DOWN_SCALING_FACTOR_X,
3763 DRIZZLE_DOWN_SCALING_FACTOR_Y,
3764 DRIZZLE_DOWN_SCALING_FACTOR_Z,
3765 aMethod, use_errorweights);
3766 hdrl_resample_outgrid_parameter *aParams_outputgrid =
3767 (hdrl_resample_outgrid_parameter *) outputgrid;
3768 hdrl_resample_util_fill_outputgrid(aParams_outputgrid);
3769 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
3770 aParams_outputgrid->wcs = hdrl_resample_util_crea_wcs_3d();
3771
3772 /* generates a cube to be resampled */
3773 hdrl_resample_result * cube = hdrl_resample_compute(ResTable, aParams_method, outputgrid,
3774 aParams_outputgrid->wcs);
3775
3776 /* test invalid input */
3777 hdrl_resample_cube_nearest(NULL, ResTable, aGrid, aParams_outputgrid);
3778 cpl_test_error(CPL_ERROR_NULL_INPUT);
3779
3780 hdrl_resample_cube_nearest(cube, NULL, aGrid, aParams_outputgrid);
3781 cpl_test_error(CPL_ERROR_NULL_INPUT);
3782
3783 hdrl_resample_cube_nearest(cube, ResTable, NULL, aParams_outputgrid);
3784 cpl_test_error(CPL_ERROR_NULL_INPUT);
3785
3786 hdrl_resample_cube_nearest(cube, ResTable, aGrid, NULL);
3787 cpl_test_error(CPL_ERROR_NULL_INPUT);
3788
3789 /* test invalid input */
3790 hdrl_resample_cube_nearest(cube, ResTable, aGrid,aParams_outputgrid);
3791 cpl_test_error(CPL_ERROR_NONE);
3792
3793 /* clean memory */
3794 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
3795 hdrl_resample_pixgrid_delete(aGrid);
3796 cpl_table_delete(ResTable);
3797 hdrl_parameter_delete(outputgrid);
3798 hdrl_parameter_delete(aParams_method);
3799 hdrl_resample_result_delete(cube);
3800 return cpl_error_get_code();
3801}
3802/*----------------------------------------------------------------------------*/
3807/*----------------------------------------------------------------------------*/
3808static cpl_error_code
3809test_hdrl_resample_cube_weighted(void) {
3810
3811 cpl_table *ResTable = NULL;
3812 const cpl_size size = 10;
3813 cpl_size x = 2;
3814 cpl_size y = 2;
3815 cpl_size z = 2;
3816 unsigned short aNMaps = 1;
3817
3818 /* prepare re-sample table */
3819 hdrl_resample_create_table(&ResTable, size);
3820
3821 /* prepare a Grid */
3822 hdrl_resample_pixgrid *aGrid = hdrl_resample_pixgrid_new(x, y, z, aNMaps);
3823
3824 /* prepare output-grid parameter */
3825 hdrl_parameter * outputgrid =
3826 hdrl_resample_parameter_create_outgrid3D_userdef(HDRL_CDELT1,
3827 HDRL_CDELT2, HDRL_CDELT3,
3828 HDRL_RA_MIN, HDRL_RA_MAX, HDRL_DEC_MIN, HDRL_DEC_MAX,
3829 HDRL_LAMBDA_MIN, HDRL_LAMBDA_MAX,5.);
3830
3831 /* creates the cube */
3832 hdrl_resample_method aMethod = HDRL_RESAMPLE_METHOD_DRIZZLE;
3833 cpl_boolean use_errorweights = CPL_TRUE;
3834 hdrl_parameter * aParams_method =
3835 hdrl_resample_util_methodparam_create (
3836 LOOP_DISTANCE,
3837 RENKA_CRITICAL_RADIUS,
3838 LANCZOS_KERNEL_SIZE,
3839 DRIZZLE_DOWN_SCALING_FACTOR_X,
3840 DRIZZLE_DOWN_SCALING_FACTOR_Y,
3841 DRIZZLE_DOWN_SCALING_FACTOR_Z,
3842 aMethod, use_errorweights);
3843 hdrl_resample_outgrid_parameter *aParams_outputgrid =
3844 (hdrl_resample_outgrid_parameter *) outputgrid;
3845 hdrl_resample_util_fill_outputgrid(aParams_outputgrid);
3846 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
3847 aParams_outputgrid->wcs = hdrl_resample_util_crea_wcs_3d();
3848
3849 hdrl_resample_method_parameter * mp =
3850 (hdrl_resample_method_parameter *) aParams_method;
3851 cpl_test_nonnull(mp);
3852 /* generates a cube to be resampled */
3853 hdrl_resample_result * cube = hdrl_resample_compute(ResTable, aParams_method, outputgrid,
3854 aParams_outputgrid->wcs);
3855
3856 /* test invalid input */
3857 hdrl_resample_cube_weighted(NULL, ResTable, aGrid, mp,
3858 aParams_outputgrid);
3859 cpl_test_error(CPL_ERROR_NULL_INPUT);
3860
3861 hdrl_resample_cube_weighted(cube, NULL, aGrid, mp,
3862 aParams_outputgrid);
3863 cpl_test_error(CPL_ERROR_NULL_INPUT);
3864
3865 hdrl_resample_cube_weighted(cube, ResTable, NULL, mp,
3866 aParams_outputgrid);
3867 cpl_test_error(CPL_ERROR_NULL_INPUT);
3868
3869 hdrl_resample_cube_weighted(cube, ResTable, aGrid, NULL,
3870 aParams_outputgrid);
3871 cpl_test_error(CPL_ERROR_NULL_INPUT);
3872
3873 hdrl_resample_cube_weighted(cube, ResTable, aGrid, mp,
3874 NULL);
3875 cpl_test_error(CPL_ERROR_NULL_INPUT);
3876
3877 /* test valid input */
3878 hdrl_resample_cube_weighted(cube, ResTable, aGrid, mp,
3879 aParams_outputgrid);
3880 cpl_test_error(CPL_ERROR_NONE);
3881
3882 mp->loop_distance = -1;
3883 mp->lanczos_kernel_size = 0;
3884 hdrl_resample_cube_weighted(cube, ResTable, aGrid, mp,
3885 aParams_outputgrid);
3886 cpl_test_error(CPL_ERROR_NONE);
3887
3888 /* clean memory */
3889 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
3890 hdrl_resample_pixgrid_delete(aGrid);
3891 cpl_table_delete(ResTable);
3892 hdrl_parameter_delete(outputgrid);
3893 hdrl_parameter_delete(aParams_method);
3894 hdrl_resample_result_delete(cube);
3895
3896 return cpl_error_get_code();
3897}
3898/*----------------------------------------------------------------------------*/
3903/*----------------------------------------------------------------------------*/
3904static cpl_error_code
3905test_hdrl_resample_cube(void) {
3906 cpl_table *ResTable = NULL;
3907 const cpl_size size = 10;
3908
3909 /* prepare re-sample table */
3910 hdrl_resample_create_table(&ResTable, size);
3911
3912 /* prepare a Grid */
3913 hdrl_resample_pixgrid *aGrid = NULL;
3914 /* prepare output-grid parameter */
3915 hdrl_parameter * outputgrid =
3916 hdrl_resample_parameter_create_outgrid3D_userdef(HDRL_CDELT1,
3917 HDRL_CDELT2, HDRL_CDELT3,
3918 HDRL_RA_MIN, HDRL_RA_MAX, HDRL_DEC_MIN, HDRL_DEC_MAX,
3919 HDRL_LAMBDA_MIN, HDRL_LAMBDA_MAX,5.);
3920
3921 /* creates the cube */
3922 hdrl_resample_method aMethod = HDRL_RESAMPLE_METHOD_DRIZZLE;
3923 cpl_boolean use_errorweights = CPL_TRUE;
3924 hdrl_parameter * aParams_method =
3925 hdrl_resample_util_methodparam_create (
3926 LOOP_DISTANCE,
3927 RENKA_CRITICAL_RADIUS,
3928 LANCZOS_KERNEL_SIZE,
3929 DRIZZLE_DOWN_SCALING_FACTOR_X,
3930 DRIZZLE_DOWN_SCALING_FACTOR_Y,
3931 DRIZZLE_DOWN_SCALING_FACTOR_Z,
3932 aMethod, use_errorweights);
3933 hdrl_resample_outgrid_parameter *aParams_outputgrid =
3934 (hdrl_resample_outgrid_parameter *) outputgrid;
3935 hdrl_resample_util_fill_outputgrid(aParams_outputgrid);
3936 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
3937 aParams_outputgrid->wcs = hdrl_resample_util_crea_wcs_3d();
3938
3939 hdrl_resample_method_parameter * mp =
3940 (hdrl_resample_method_parameter *) aParams_method;
3941 cpl_test_nonnull(mp);
3942 /* test invalid input */
3943 hdrl_resample_cube(NULL, mp, aParams_outputgrid, &aGrid);
3944 cpl_test_error(CPL_ERROR_NULL_INPUT);
3945
3946 hdrl_resample_cube(ResTable, NULL, aParams_outputgrid, &aGrid);
3947 cpl_test_error(CPL_ERROR_NULL_INPUT);
3948
3949 hdrl_resample_cube(ResTable, mp, NULL, &aGrid);
3950 cpl_test_error(CPL_ERROR_NULL_INPUT);
3951
3952 hdrl_resample_cube(ResTable, mp, aParams_outputgrid, NULL);
3953 cpl_test_error(CPL_ERROR_NULL_INPUT);
3954
3955 /* test valid input */
3956 hdrl_resample_result * cube = hdrl_resample_cube(ResTable, mp, aParams_outputgrid, &aGrid);
3957 cpl_test_error(CPL_ERROR_NONE);
3958 hdrl_resample_result_delete(cube);
3959 hdrl_resample_pixgrid_delete(aGrid);
3960
3961 /* test wrong method input */
3962 mp->method = 10;
3963 cube = hdrl_resample_cube(ResTable, mp, aParams_outputgrid, &aGrid);
3964 cpl_test_error(CPL_ERROR_UNSUPPORTED_MODE);
3965 hdrl_resample_result_delete(cube);
3966 cpl_table_delete(ResTable);
3967 hdrl_resample_pixgrid_delete(aGrid);
3968
3969 /* test invalid grid creation input */
3970 ResTable = cpl_table_new(0);
3971 cube = hdrl_resample_cube(ResTable, mp, aParams_outputgrid, &aGrid);
3972 cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
3973
3974 /* clean memory */
3975 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
3976 hdrl_resample_pixgrid_delete(aGrid);
3977 cpl_table_delete(ResTable);
3978 hdrl_parameter_delete(outputgrid);
3979 hdrl_parameter_delete(aParams_method);
3980 hdrl_resample_result_delete(cube);
3981
3982 cpl_test_error(CPL_ERROR_NONE);
3983 return cpl_error_get_code();
3984}
3985/*----------------------------------------------------------------------------*/
3990/*----------------------------------------------------------------------------*/
3991static cpl_error_code
3992test_hdrl_resampling_set_outputgrid(void)
3993{
3994 int xsize = HDRL_SIZE_X;
3995 int ysize = HDRL_SIZE_Y;
3996 int zsize = HDRL_SIZE_Z;
3997
3998 cpl_table *ResTable = NULL;
3999 const cpl_size size = 10;
4000
4001 /* prepare re-sample table */
4002 hdrl_resample_create_table(&ResTable, size);
4003
4004 /* prepare a Grid */
4005 hdrl_resample_pixgrid *aGrid = NULL;
4006
4007 /* prepare output-grid parameter */
4008 hdrl_parameter * outputgrid =
4009 hdrl_resample_parameter_create_outgrid3D_userdef(HDRL_CDELT1,
4010 HDRL_CDELT2, HDRL_CDELT3,
4011 HDRL_RA_MIN, HDRL_RA_MAX, HDRL_DEC_MIN, HDRL_DEC_MAX,
4012 HDRL_LAMBDA_MIN, HDRL_LAMBDA_MAX,5.);
4013
4014 /* creates the cube */
4015 hdrl_resample_method aMethod = HDRL_RESAMPLE_METHOD_DRIZZLE;
4016 cpl_boolean use_errorweights = CPL_TRUE;
4017 hdrl_parameter * aParams_method =
4018 hdrl_resample_util_methodparam_create (
4019 LOOP_DISTANCE,
4020 RENKA_CRITICAL_RADIUS,
4021 LANCZOS_KERNEL_SIZE,
4022 DRIZZLE_DOWN_SCALING_FACTOR_X,
4023 DRIZZLE_DOWN_SCALING_FACTOR_Y,
4024 DRIZZLE_DOWN_SCALING_FACTOR_Z,
4025 aMethod, use_errorweights);
4026 hdrl_resample_outgrid_parameter *aParams_outputgrid =
4027 (hdrl_resample_outgrid_parameter *) outputgrid;
4028 hdrl_resample_util_fill_outputgrid(aParams_outputgrid);
4029 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
4030 aParams_outputgrid->wcs = hdrl_resample_util_crea_wcs_3d();
4031
4032 hdrl_resample_method_parameter * mp =
4033 (hdrl_resample_method_parameter *) aParams_method;
4034 cpl_test_nonnull(mp);
4035
4036 hdrl_resample_result * cube = hdrl_resample_cube(ResTable, mp,
4037 aParams_outputgrid, &aGrid);
4038
4039 /* test invalid input */
4040 hdrl_resampling_set_outputgrid (-1, ysize, zsize, cube, aParams_outputgrid);
4041 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
4042
4043 hdrl_resampling_set_outputgrid (xsize, -1, zsize, cube, aParams_outputgrid);
4044 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
4045
4046 hdrl_resampling_set_outputgrid (xsize, ysize, -1, cube, aParams_outputgrid);
4047 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
4048
4049 hdrl_resampling_set_outputgrid (xsize, ysize, zsize, NULL, aParams_outputgrid);
4050 cpl_test_error(CPL_ERROR_NULL_INPUT);
4051
4052 hdrl_resampling_set_outputgrid (xsize, ysize, zsize, cube, NULL);
4053 cpl_test_error(CPL_ERROR_NULL_INPUT);
4054 /* test valid input */
4055 /* We need to do delete the header cube->header, as
4056 * hdrl_resampling_set_outputgrid allocates mem for cube->header
4057 * even if cube->header might already exists
4058 */
4059 cpl_propertylist_delete(cube->header);
4060
4061 hdrl_resampling_set_outputgrid (xsize, ysize, zsize, cube, aParams_outputgrid);
4062 cpl_test_error(CPL_ERROR_NONE);
4063
4064 cpl_wcs_delete((cpl_wcs *)aParams_outputgrid->wcs);
4065 hdrl_resample_pixgrid_delete(aGrid);
4066 cpl_table_delete(ResTable);
4067 hdrl_parameter_delete(outputgrid);
4068 hdrl_parameter_delete(aParams_method);
4069 hdrl_resample_result_delete(cube);
4070 cpl_test_error(CPL_ERROR_NONE);
4071 return cpl_error_get_code();
4072}
4073static cpl_error_code
4074test_hdrl_resample_inputtable_verify(void)
4075{
4076 cpl_table* tab = NULL;
4077 hdrl_resample_inputtable_verify(NULL);
4078 cpl_test_error(CPL_ERROR_NULL_INPUT);
4079
4080 /* verify column existence */
4081 tab= cpl_table_new(1);
4082 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_RA, CPL_TYPE_DOUBLE);
4083 hdrl_resample_inputtable_verify(tab);
4084 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4085
4086 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_DEC, CPL_TYPE_DOUBLE);
4087 hdrl_resample_inputtable_verify(tab);
4088 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4089
4090 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_LAMBDA, CPL_TYPE_DOUBLE);
4091 hdrl_resample_inputtable_verify(tab);
4092 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4093
4094 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_BPM, CPL_TYPE_INT);
4095 hdrl_resample_inputtable_verify(tab);
4096 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4097
4098 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_ERRORS, CPL_TYPE_DOUBLE);
4099 hdrl_resample_inputtable_verify(tab);
4100 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4101
4102 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_DATA, CPL_TYPE_DOUBLE);
4103 hdrl_resample_inputtable_verify(tab);
4104 cpl_test_error(CPL_ERROR_NONE);
4105 cpl_table_delete(tab);
4106
4107 /* verify column type */
4108 tab= cpl_table_new(1);
4109 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_RA, CPL_TYPE_INT);
4110 hdrl_resample_inputtable_verify(tab);
4111 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4112 cpl_table_erase_column(tab, HDRL_RESAMPLE_TABLE_RA);
4113 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_RA, CPL_TYPE_DOUBLE);
4114
4115 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_DEC, CPL_TYPE_INT);
4116 hdrl_resample_inputtable_verify(tab);
4117 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4118 cpl_table_erase_column(tab, HDRL_RESAMPLE_TABLE_DEC);
4119 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_DEC, CPL_TYPE_DOUBLE);
4120
4121 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_LAMBDA, CPL_TYPE_INT);
4122 hdrl_resample_inputtable_verify(tab);
4123 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4124 cpl_table_erase_column(tab, HDRL_RESAMPLE_TABLE_LAMBDA);
4125 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_LAMBDA, CPL_TYPE_DOUBLE);
4126
4127 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_BPM, CPL_TYPE_DOUBLE);
4128 hdrl_resample_inputtable_verify(tab);
4129 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4130 cpl_table_erase_column(tab, HDRL_RESAMPLE_TABLE_BPM);
4131 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_BPM, CPL_TYPE_INT);
4132
4133 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_ERRORS, CPL_TYPE_INT);
4134 hdrl_resample_inputtable_verify(tab);
4135 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4136 cpl_table_erase_column(tab, HDRL_RESAMPLE_TABLE_ERRORS);
4137 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_ERRORS, CPL_TYPE_DOUBLE);
4138
4139 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_DATA, CPL_TYPE_INT);
4140 hdrl_resample_inputtable_verify(tab);
4141 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
4142 cpl_table_erase_column(tab, HDRL_RESAMPLE_TABLE_DATA);
4143 cpl_table_new_column(tab, HDRL_RESAMPLE_TABLE_DATA, CPL_TYPE_DOUBLE);
4144 hdrl_resample_inputtable_verify(tab);
4145 cpl_test_error(CPL_ERROR_NONE);
4146
4147 cpl_table_delete(tab);
4148 return cpl_error_get_code();
4149}
4150
4154static cpl_propertylist*
4155hdrl_resample_util_header_create(const int sx, const int sy,
4156 const double dx, const double dy)
4157{
4158
4159 cpl_propertylist* plist = cpl_propertylist_new();
4160 cpl_propertylist_append_bool(plist,"SIMPLE", CPL_TRUE);
4161 cpl_propertylist_append_int(plist,"NAXIS", 2);
4162 cpl_propertylist_append_int(plist,"NAXIS1", sx);
4163 cpl_propertylist_append_int(plist,"NAXIS2", sy);
4164 cpl_propertylist_append_bool(plist,"EXTEND", CPL_TRUE);
4165 cpl_propertylist_append_double(plist,"CRPIX1", 0.5 * sx);
4166 cpl_propertylist_append_double(plist,"CRPIX2", 0.5 * sy);
4167 cpl_propertylist_append_double(plist,"CRVAL1", 0.5 * sx);
4168 cpl_propertylist_append_double(plist,"CRVAL2", 0.5 * sy);
4169 cpl_propertylist_append_double(plist,"CDELT1", dx );
4170 cpl_propertylist_append_double(plist,"CDELT2", dy );
4171 cpl_propertylist_append_double(plist,"CD1_1", -dx );
4172 cpl_propertylist_append_double(plist,"CD1_2", 0. );
4173 cpl_propertylist_append_double(plist,"CD2_1", 0. );
4174 cpl_propertylist_append_double(plist,"CD2_2", dy );
4175 cpl_propertylist_append_string(plist,"CTYPE1","PIXEL");
4176 cpl_propertylist_append_string(plist,"CTYPE2","PIXEL");
4177 cpl_propertylist_append_string(plist,"CUNIT1","PIXEL");
4178 cpl_propertylist_append_string(plist,"CUNIT2","PIXEL");
4179
4180 return plist;
4181}
4182
4186static cpl_image*
4187hdrl_resample_util_crea_image_uniform(const cpl_size sx, const cpl_size sy,
4188 const double value)
4189{
4190 assert(value > 0);
4191 cpl_image * image = cpl_image_new(sx, sy, CPL_TYPE_DOUBLE);
4192 cpl_image_add_scalar(image,value);
4193 return image;
4194}
4195
4196static double
4197hdrldemo_get_resampled_pix_value(const char* type, const double pix_value,
4198 const double outlier, cpl_boolean is_bad,
4199 hdrl_resample_method aMethod)
4200{
4201 cpl_size sx = HDRL_SIZE_X;
4202 cpl_size sy = HDRL_SIZE_X;
4203 cpl_size xc = 0.5 * sx;
4204 cpl_size yc = 0.5 * sy;
4205
4206 const double value = pix_value;
4207 const double dx = 0.01;
4208 const double dy = 0.01;
4209 cpl_image* simul = NULL;
4210 cpl_image* errs = NULL;
4211
4212 cpl_test_error(CPL_ERROR_NONE);
4213 if(strcmp(type,"uniform") == 0) {
4214 simul = hdrl_resample_util_crea_image_uniform(sx, sy, value);
4215 } else {
4216 simul = hdrl_resample_util_crea_image_dice_5(sx, sy, value);
4217 }
4218 cpl_test_error(CPL_ERROR_NONE);
4219
4220 char* fname;
4221 if( is_bad ) {
4222 cpl_image_set(simul, xc, yc, outlier);
4223 fname = cpl_sprintf("cube_and_bp.fits");
4224 } else {
4225 fname = cpl_sprintf("cube_not_bp.fits");
4226 }
4227 cpl_free(fname);
4228 cpl_test_error(CPL_ERROR_NONE);
4229 errs = cpl_image_duplicate(simul);
4230 cpl_image_power(errs, 0.5);
4231
4232 /* We need to add wcs */
4233 cpl_propertylist* plist = NULL;
4234 plist = hdrl_resample_util_header_create(sx, sy, dx, dy);
4235 cpl_image_save(simul, "image.fits", CPL_TYPE_DOUBLE,plist,CPL_IO_DEFAULT);
4236 cpl_propertylist_delete(plist);
4237 plist = cpl_propertylist_load("image.fits",0);
4238 cpl_wcs *wcs = cpl_wcs_new_from_propertylist(plist);
4239 cpl_test_nonnull(wcs);
4240 cpl_test_error(CPL_ERROR_NONE);
4241
4242 /* we now create the MUSE table */
4243 cpl_imagelist* ilist = NULL;
4244 cpl_imagelist* elist = NULL;
4245 ilist = cpl_imagelist_new();
4246 elist = cpl_imagelist_new();
4247
4248 cpl_imagelist_set(ilist, simul, 0);
4249 cpl_imagelist_set(elist, errs, 0);
4250 cpl_table* pixel_table = NULL;
4251 hdrl_imagelist* hlist = NULL;
4252 if( is_bad == CPL_TRUE) {
4253 cpl_image* qual = cpl_image_new(sx, sy, CPL_TYPE_INT);
4254 cpl_image_set(qual, xc, yc, 1);
4255 cpl_imagelist* qlist = cpl_imagelist_new();
4256 cpl_imagelist_set(qlist, qual, 0);
4257 hlist = hdrl_imagelist_create(ilist, elist);
4258 pixel_table = hdrl_resample_imagelist_to_table(hlist, wcs);
4259 cpl_imagelist_delete(qlist);
4260 } else {
4261 hlist = hdrl_imagelist_create(ilist, elist);
4262 pixel_table = hdrl_resample_imagelist_to_table(hlist, wcs);
4263 }
4264 cpl_imagelist_delete(ilist);
4265 cpl_imagelist_delete(elist);
4266 hdrl_imagelist_delete(hlist);
4267
4268 cpl_frameset* frameset = cpl_frameset_new();
4269
4270 /* now we resample */
4271 cpl_msg_info(cpl_func,"start resample");
4272
4273 cpl_test_error(CPL_ERROR_NONE);
4274 cpl_wcs_delete(wcs);
4275 cpl_propertylist* p = cpl_propertylist_new();
4276 cpl_propertylist_append_string(p,"CTYPE1", "RA---TAN");
4277 cpl_propertylist_append_string(p,"CTYPE2", "DEC--TAN");
4278 cpl_propertylist_append_double(p,"CRVAL1", 0.);
4279 cpl_propertylist_append_double(p,"CRVAL2", 0.);
4280 cpl_propertylist_append_double(p,"CRPIX1", xc);
4281 cpl_propertylist_append_double(p,"CRPIX2", yc);
4282 cpl_propertylist_append_double(p,"CD1_1", -dx);
4283 cpl_propertylist_append_double(p,"CD1_2", 0.1); /* add rotation if not 0 */
4284 cpl_propertylist_append_double(p,"CD2_1", 0.1); /* add rotation if not 0 */
4285 cpl_propertylist_append_double(p,"CD2_2", dy);
4286 wcs = cpl_wcs_new_from_propertylist(p);
4287 cpl_test_nonnull(wcs);
4288 cpl_propertylist_delete(p);
4289
4290 const double ramin = xc - 0.5*sx*dx;
4291 const double ramax = xc + 0.5*sx*dx;
4292 const double decmin = yc - 0.5*sy*dy;
4293 const double decmax = yc + 0.5*sy*dy;
4294 const double lambmin = 550;
4295 const double lambmax = 551;
4296 const double dlambda = 1;
4297
4298 hdrl_parameter* aParams_outputgrid =
4299 hdrl_resample_parameter_create_outgrid3D_userdef(dx, dy, dlambda,
4300 ramin, ramax, decmin, decmax, lambmin, lambmax, 5.);
4301
4302 cpl_boolean use_errorweights = CPL_FALSE;
4303 hdrl_parameter * aParams_method =
4304 hdrl_resample_util_methodparam_create (
4305 LOOP_DISTANCE,
4306 RENKA_CRITICAL_RADIUS,
4307 LANCZOS_KERNEL_SIZE,
4308 DRIZZLE_DOWN_SCALING_FACTOR_X,
4309 DRIZZLE_DOWN_SCALING_FACTOR_Y,
4310 DRIZZLE_DOWN_SCALING_FACTOR_Z,
4311 aMethod, use_errorweights);
4312
4313 hdrl_resample_result *cube = NULL;
4314 cube = hdrl_resample_compute(pixel_table, aParams_method, aParams_outputgrid, wcs);
4315 cpl_test_error(CPL_ERROR_NONE);
4316
4317 hdrl_imagelist* himlist = cube->himlist;
4318 hdrl_image* hima = hdrl_imagelist_get(himlist, 0);
4319 double* data = cpl_image_get_data(hdrl_image_get_image(hima));
4320 cpl_size rsx = hdrl_imagelist_get_size_x(himlist);
4321 cpl_size rsy = hdrl_imagelist_get_size_y(himlist);
4322 cpl_size rpix = 0.5 * rsy * rsx + rsx;
4323
4324 double pix_ref = data[rpix];
4325 cpl_test_error(CPL_ERROR_NONE);
4326 cpl_test_error(CPL_ERROR_NONE);
4327
4328 /* clean memory */
4329 cpl_free(aParams_method);
4330 cpl_wcs_delete(wcs);
4331 hdrl_parameter_destroy(aParams_outputgrid);
4332 cpl_table_delete(pixel_table);
4333 hdrl_resample_result_delete(cube);
4334 cpl_frameset_delete(frameset);
4335 cpl_propertylist_delete(plist);
4336 return pix_ref;
4337}
4338
4339static cpl_error_code
4340test_resample_image_dice(hdrl_resample_method aMethod)
4341{
4342
4343 /* Create an image with the number 5 of a dice.
4344 Central point is round
4345 Other points are elliptical with major axis oriented along different
4346 directions
4347 */
4348 cpl_test_error(CPL_ERROR_NONE);
4349 double ref_value =
4350 hdrldemo_get_resampled_pix_value("uniform", HDRL_FLUX_ADU, 1000., CPL_FALSE,
4351 aMethod );
4352
4353 double check_value =
4354 hdrldemo_get_resampled_pix_value("uniform", HDRL_FLUX_ADU, 100000., CPL_TRUE,
4355 aMethod);
4356
4357 cpl_test_abs(check_value, ref_value, HDRL_EPS_TEST);
4358 /* possible further tests:
4359 move the image horizontally
4360 move the image vertically
4361 move the image diagonally
4362 rotate the image around a point: around the center, at the edge
4363 over-sample the image
4364 sub-sample the image
4365 */
4366 return cpl_error_get_code();
4367}
4368
4369static hdrl_resample_result *
4370hdrl_resample_util_get_resampled_pix_cube(const char* type, const double pix_value,
4371 const double outlier, cpl_boolean is_bad, const char* suffix)
4372{
4373 cpl_size sx = HDRL_SIZE_X;
4374 cpl_size sy = HDRL_SIZE_Y;
4375 cpl_size xc = 0.5 * sx;
4376 cpl_size yc = 0.5 * sy;
4377 const double value = pix_value;
4378 const double dx = 0.01;
4379 const double dy = 0.01;
4380 cpl_image* simul = NULL;
4381 cpl_image* errs = NULL;
4382
4383 cpl_test_error(CPL_ERROR_NONE);
4384 if(strcmp(type,"uniform") == 0) {
4385 simul = hdrl_resample_util_crea_image_uniform(sx, sy, value);
4386 } else {
4387 simul = hdrl_resample_util_crea_image_dice_5(sx, sy, value);
4388 }
4389 cpl_test_error(CPL_ERROR_NONE);
4390
4391 char* bname;
4392 if( is_bad ) {
4393 cpl_msg_warning(cpl_func,"is bad");
4394 cpl_image_set(simul, xc, yc, outlier);
4395 bname = cpl_sprintf("%s_%s.fits","cube_and_bp",suffix);
4396 } else {
4397 cpl_msg_warning(cpl_func,"is not bad");
4398 bname = cpl_sprintf("%s_%s.fits","cube_not_bp",suffix);
4399 }
4400 cpl_test_error(CPL_ERROR_NONE);
4401
4402 errs = cpl_image_duplicate(simul);
4403 cpl_image_power(errs, 0.5);
4404
4405 /* Add wcs */
4406 cpl_propertylist* plist = NULL;
4407 plist = hdrl_resample_util_header_create(sx, sy, dx, dy);
4408 cpl_image_save(simul, "image.fits", CPL_TYPE_DOUBLE,plist,CPL_IO_DEFAULT);
4409 cpl_propertylist_delete(plist);
4410 plist = cpl_propertylist_load("image.fits",0);
4411 cpl_wcs *wcs = cpl_wcs_new_from_propertylist(plist);
4412 cpl_test_nonnull(wcs);
4413
4414 /* we now create the pixel table */
4415 cpl_imagelist* ilist = NULL;
4416 cpl_imagelist* elist = NULL;
4417 cpl_table* pixel_table = NULL;
4418 ilist = cpl_imagelist_new();
4419 elist = cpl_imagelist_new();
4420 cpl_imagelist_set(ilist, simul, 0);
4421 cpl_imagelist_set(elist, errs, 0);
4422 cpl_test_error(CPL_ERROR_NONE);
4423 cpl_imagelist_save(ilist,bname, CPL_TYPE_DOUBLE, NULL, CPL_IO_DEFAULT);
4424 cpl_imagelist_save(elist,bname, CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
4425 cpl_test_error(CPL_ERROR_NONE);
4426 hdrl_imagelist* hlist;
4427 if( is_bad == CPL_TRUE) {
4428 cpl_image* qual = cpl_image_new(sx, sy, CPL_TYPE_INT);
4429 cpl_image_set(qual, xc, yc, 1);
4430 cpl_imagelist* qlist = cpl_imagelist_new();
4431 cpl_imagelist_set(qlist, qual, 0);
4432 cpl_imagelist_save(qlist,bname, CPL_TYPE_INT, NULL, CPL_IO_EXTEND);
4433 hlist = hdrl_imagelist_create(ilist, elist);
4434 pixel_table = hdrl_resample_imagelist_to_table(hlist, wcs);
4435 cpl_imagelist_delete(qlist);
4436 } else {
4437 hlist = hdrl_imagelist_create(ilist, elist);
4438 pixel_table = hdrl_resample_imagelist_to_table(hlist, wcs);
4439 }
4440 cpl_test_error(CPL_ERROR_NONE);
4441
4442 cpl_imagelist_delete(ilist);
4443 cpl_imagelist_delete(elist);
4444 hdrl_imagelist_delete(hlist);
4445 cpl_free(bname);
4446 cpl_frameset* frameset = cpl_frameset_new();
4447
4448 /* now we resample */
4449 cpl_msg_info(cpl_func,"start resample");
4450 const double ramin = xc - 0.5*sx*dx;
4451 const double ramax = xc + 0.5*sx*dx;
4452 const double decmin = yc - 0.5*sy*dy;
4453 const double decmax = yc + 0.5*sy*dy;
4454 const double lambmin = 550;
4455 const double lambmax = 551;
4456 const double dlambda = 1;
4457 cpl_wcs_delete(wcs);
4458 cpl_propertylist* p = cpl_propertylist_new();
4459 cpl_propertylist_append_string(p,"CTYPE1", "RA---TAN");
4460 cpl_propertylist_append_string(p,"CTYPE2", "DEC--TAN");
4461 cpl_propertylist_append_double(p,"CRVAL1", 0.);
4462 cpl_propertylist_append_double(p,"CRVAL2", 0.);
4463 cpl_propertylist_append_double(p,"CRPIX1", xc);
4464 cpl_propertylist_append_double(p,"CRPIX2", yc);
4465 cpl_propertylist_append_double(p,"CD1_1", -dx);
4466 cpl_propertylist_append_double(p,"CD1_2", 0.);
4467 cpl_propertylist_append_double(p,"CD2_1", 0.);
4468 cpl_propertylist_append_double(p,"CD2_2", dy);
4469 wcs = cpl_wcs_new_from_propertylist(p);
4470 cpl_propertylist_delete(p);
4471
4472 hdrl_parameter* aParams_outputgrid =
4473 hdrl_resample_parameter_create_outgrid3D_userdef(dx, dy, dlambda,
4474 ramin, ramax, decmin, decmax, lambmin, lambmax, 5.);
4475
4476 cpl_boolean use_errorweights = CPL_FALSE;
4477 /* RENKA, LANCZOS, DRIZZLE generate a cube with all flagged pixels */
4478 hdrl_resample_method aMethod = HDRL_RESAMPLE_METHOD_QUADRATIC;
4479 hdrl_parameter * aParams_method =
4480 hdrl_resample_util_methodparam_create (
4481 LOOP_DISTANCE,
4482 RENKA_CRITICAL_RADIUS,
4483 LANCZOS_KERNEL_SIZE,
4484 DRIZZLE_DOWN_SCALING_FACTOR_X,
4485 DRIZZLE_DOWN_SCALING_FACTOR_Y,
4486 DRIZZLE_DOWN_SCALING_FACTOR_Z,
4487 aMethod, use_errorweights);
4488
4489 hdrl_resample_result *cube = NULL;
4490 cube = hdrl_resample_compute(pixel_table, aParams_method,
4491 aParams_outputgrid, wcs);
4492
4493 cpl_msg_info(cpl_func,"end resample");
4494 /* clean memory */
4495 cpl_table_delete(pixel_table);
4496 cpl_wcs_delete(wcs);
4497 hdrl_parameter_destroy(aParams_outputgrid);
4498 hdrl_parameter_destroy(aParams_method);
4499 cpl_frameset_delete(frameset);
4500 cpl_propertylist_delete(plist);
4501
4502 return cube;
4503}
4504
4505static cpl_error_code
4506test_resample_image_with_outlier(void)
4507{
4508 cpl_test_error(CPL_ERROR_NONE);
4509
4510 hdrl_resample_result *cube1 =
4511 hdrl_resample_util_get_resampled_pix_cube("uniform", HDRL_FLUX_ADU, 1000., CPL_FALSE,"one");
4512 hdrl_image* hima1 = hdrl_imagelist_get(cube1->himlist,0);
4513
4514 hdrl_resample_result *cube2 =
4515 hdrl_resample_util_get_resampled_pix_cube("uniform", HDRL_FLUX_ADU, 100000., CPL_FALSE,"two");
4516
4517 hdrl_image* hima2 = hdrl_imagelist_get(cube2->himlist,0);
4518 hdrl_image* hdiff = hdrl_image_sub_image_create(hima1, hima2);
4519 cpl_image* diff = hdrl_image_get_image(hdiff);
4520
4521 cpl_image_save(diff,"diff.fits",CPL_TYPE_DOUBLE,NULL,CPL_IO_DEFAULT);
4522 cpl_test_error(CPL_ERROR_NONE);
4523
4524 cpl_test_abs(0,cpl_image_get_mean(diff), HDRL_EPS_TEST);
4525 cpl_test_abs(0,cpl_image_get_stdev(diff), HDRL_EPS_TEST);
4526 cpl_test_abs(0,cpl_image_get_min(diff), HDRL_EPS_TEST);
4527 cpl_test_abs(0,cpl_image_get_max(diff), HDRL_EPS_TEST);
4528
4529 hdrl_image_delete(hdiff);
4530 hdrl_resample_result_delete(cube1);
4531 hdrl_resample_result_delete(cube2);
4532 cpl_test_error(CPL_ERROR_NONE);
4533 return cpl_error_get_code();
4534}
4535
4536/*----------------------------------------------------------------------------*/
4540/*----------------------------------------------------------------------------*/
4541int main(void)
4542{
4543 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);
4544
4545 test_resample_wcs_as_muse();
4546 test_hdrl_resample_wcs_print();
4547 test_hdrl_resample_outgrid_parameter_print();
4548 test_hdrl_resample_method_parameter_print();
4549 test_hdrl_wcs_xy_to_radec();
4550 test_hdrl_resample_pfits_get();
4551 test_hdrl_resample_smallwcs_new();
4552 test_hdrl_resample_pixgrid_delete();
4553 test_hdrl_resample_compute();
4554
4555 test_resample_parameters();
4556 test_hdrl_wcs_to_propertylist();
4557 test_hdrl_resample_image_to_table();
4558 test_hdrl_resample_imagelist_to_table();
4559 test_hdrl_resample_parameter_create_outgrid2D();
4560 test_hdrl_resample_parameter_create_outgrid3D();
4561 test_hdrl_resample_parameter_create_outgrid2D_userdef();
4562 test_hdrl_resample_parameter_create_outgrid3D_userdef();
4563
4564 test_hdrl_resample_weight_function_renka();
4565 test_hdrl_resample_weight_function_linear();
4566 test_hdrl_resample_weight_function_quadratic();
4567 test_hdrl_resample_weight_function_sinc();
4568 test_hdrl_resample_weight_function_lanczos();
4569 test_hdrl_resample_weight_function_drizzle();
4570
4571 test_hdrl_resample_wcs_projplane_from_celestial();
4572 test_hdrl_resample_wcs_pixel_from_celestial_fast();
4573 test_hdrl_resample_compute_size();
4574 test_hdrl_resample_compute_size_rotated_wcs_bug();
4575 test_hdrl_resample_pixgrid_add();
4576 test_hdrl_resample_pixgrid_get_count();
4577 test_hdrl_resample_pixgrid_get_index();
4578 test_hdrl_resample_pixgrid_new();
4579
4580 test_hdrl_resample_wcs_get_scales();
4581 test_hdrl_resample_create_table();
4582 test_hdrl_resample_pixgrid_get_rows();
4583 test_hdrl_resample_pixgrid_create();
4584 test_hdrl_resample_cube_nearest();
4585 test_hdrl_resample_cube_weighted();
4586 test_hdrl_resample_cube();
4587
4588 test_hdrl_resampling_set_outputgrid();
4589 test_hdrl_resample_compute2D_multiple();
4590 test_hdrl_resample_compute3D_multiple();
4591 test_hdrl_resample_inputtable_verify();
4592
4593 test_resample_image_dice(HDRL_RESAMPLE_METHOD_LINEAR);
4594 test_resample_image_dice(HDRL_RESAMPLE_METHOD_QUADRATIC);
4595 test_resample_image_dice(HDRL_RESAMPLE_METHOD_NEAREST);
4596 test_resample_image_dice(HDRL_RESAMPLE_METHOD_RENKA);
4597 test_resample_image_dice(HDRL_RESAMPLE_METHOD_DRIZZLE);
4598 test_resample_image_dice(HDRL_RESAMPLE_METHOD_LANCZOS);
4599
4600 test_resample_image_with_outlier();
4601 return cpl_test_end(0);
4602}
4603
hdrl_image_sub_image_create
hdrl_image * hdrl_image_sub_image_create(const hdrl_image *self, const hdrl_image *other)
Subtract two images.
Definition hdrl_image_math.c:158
hdrl_image_get_mask
cpl_mask * hdrl_image_get_mask(hdrl_image *himg)
get cpl bad pixel mask from image
Definition hdrl_image.c:157
hdrl_image_get_error
cpl_image * hdrl_image_get_error(hdrl_image *himg)
get error as cpl image
Definition hdrl_image.c:131
hdrl_image_get_mask_const
const cpl_mask * hdrl_image_get_mask_const(const hdrl_image *himg)
get cpl bad pixel mask from image
Definition hdrl_image.c:175
hdrl_image_get_error_const
const cpl_image * hdrl_image_get_error_const(const hdrl_image *himg)
get error as cpl image
Definition hdrl_image.c:144
hdrl_image_create
hdrl_image * hdrl_image_create(const cpl_image *image, const cpl_image *error)
create a new hdrl_image from to existing images by copying them
Definition hdrl_image.c:295
hdrl_image_get_image
cpl_image * hdrl_image_get_image(hdrl_image *himg)
get data as cpl image
Definition hdrl_image.c:105
hdrl_image_get_image_const
const cpl_image * hdrl_image_get_image_const(const hdrl_image *himg)
get data as cpl image
Definition hdrl_image.c:118
hdrl_image_delete
void hdrl_image_delete(hdrl_image *himg)
delete hdrl_image
Definition hdrl_image.c:379
hdrl_imagelist_get_size_y
cpl_size hdrl_imagelist_get_size_y(const hdrl_imagelist *himlist)
Get number of rows of images in the imagelist.
Definition hdrl_imagelist_io.c:188
hdrl_imagelist_create
hdrl_imagelist * hdrl_imagelist_create(cpl_imagelist *imlist, cpl_imagelist *errlist)
Create an hdrl_imagelist out of 2 cpl_imagelist.
Definition hdrl_imagelist_io.c:102
hdrl_imagelist_delete
void hdrl_imagelist_delete(hdrl_imagelist *himlist)
Free all memory used by a hdrl_imagelist object including the images.
Definition hdrl_imagelist_io.c:384
hdrl_imagelist_get_const
const hdrl_image * hdrl_imagelist_get_const(const hdrl_imagelist *himlist, cpl_size inum)
Get an image from a list of images.
Definition hdrl_imagelist_io.c:229
hdrl_imagelist_get_size
cpl_size hdrl_imagelist_get_size(const hdrl_imagelist *himlist)
Get the number of images in the imagelist.
Definition hdrl_imagelist_io.c:148
hdrl_imagelist_get_size_x
cpl_size hdrl_imagelist_get_size_x(const hdrl_imagelist *himlist)
Get number of colums of images in the imagelist.
Definition hdrl_imagelist_io.c:168
hdrl_imagelist_get
hdrl_image * hdrl_imagelist_get(const hdrl_imagelist *himlist, cpl_size inum)
Get an image from a list of images.
Definition hdrl_imagelist_io.c:210
hdrl_parameter_destroy
void hdrl_parameter_destroy(hdrl_parameter *obj)
deep delete of a parameter
Definition hdrl_parameter.c:164
hdrl_parameter_delete
void hdrl_parameter_delete(hdrl_parameter *obj)
shallow delete of a parameter
Definition hdrl_parameter.c:148
hdrl_resample_parameter_create_nearest
hdrl_parameter * hdrl_resample_parameter_create_nearest(void)
Creates a resample nearest neighbor hdrl parameter object. The algorithm does not use any weighting f...
Definition hdrl_resample.c:3150
hdrl_resample_parameter_method_verify
cpl_error_code hdrl_resample_parameter_method_verify(const hdrl_parameter *hp)
verify parameters have proper values
Definition hdrl_resample.c:3381
hdrl_resample_parameter_method_check
int hdrl_resample_parameter_method_check(const hdrl_parameter *self)
check method is of proper type
Definition hdrl_resample.c:3482
hdrl_resample_parameter_create_renka
hdrl_parameter * hdrl_resample_parameter_create_renka(const int loop_distance, cpl_boolean use_errorweights, const double critical_radius)
Creates a resample renka hdrl parameter object. The algorithm uses a modified Shepard-like distance w...
Definition hdrl_resample.c:3013
hdrl_resample_parameter_create_lanczos
hdrl_parameter * hdrl_resample_parameter_create_lanczos(const int loop_distance, cpl_boolean use_errorweights, const int kernel_size)
Creates a resample Lanczos hdrl parameter object. The algorithm uses a restricted SINC distance weigh...
Definition hdrl_resample.c:3198
hdrl_resample_parameter_create_drizzle
hdrl_parameter * hdrl_resample_parameter_create_drizzle(const int loop_distance, cpl_boolean use_errorweights, const double pix_frac_x, const double pix_frac_y, const double pix_frac_lambda)
Creates a resample drizzle hdrl parameter object. The algorithm uses a drizzle-like distance weightin...
Definition hdrl_resample.c:3253
hdrl_resample_parameter_create_quadratic
hdrl_parameter * hdrl_resample_parameter_create_quadratic(const int loop_distance, cpl_boolean use_errorweights)
Creates a resample quadratic hdrl parameter object. The algorithm uses a quadratic inverse distance w...
Definition hdrl_resample.c:3108
hdrl_resample_parameter_create_linear
hdrl_parameter * hdrl_resample_parameter_create_linear(const int loop_distance, cpl_boolean use_errorweights)
Creates a resample linear hdrl parameter object. The algorithm uses a linear inverse distance weighti...
Definition hdrl_resample.c:3060
hdrl_resample_parameter_create_outgrid3D
hdrl_parameter * hdrl_resample_parameter_create_outgrid3D(const double delta_ra, const double delta_dec, const double delta_lambda)
Creates a resample_outgrid hdrl parameter object for a 3 dimensional interpolation,...
Definition hdrl_resample.c:2819
hdrl_resample_parameter_create_outgrid3D_userdef
hdrl_parameter * hdrl_resample_parameter_create_outgrid3D_userdef(const double delta_ra, const double delta_dec, const double delta_lambda, const double ra_min, const double ra_max, const double dec_min, const double dec_max, const double lambda_min, const double lambda_max, const double fieldmargin)
Creates a resample_outgrid hdrl parameter object for a 3 dimensional interpolation,...
Definition hdrl_resample.c:2948
hdrl_resample_parameter_outgrid_check
int hdrl_resample_parameter_outgrid_check(const hdrl_parameter *self)
check method is of proper type
Definition hdrl_resample.c:3465
hdrl_resample_parameter_create_outgrid2D
hdrl_parameter * hdrl_resample_parameter_create_outgrid2D(const double delta_ra, const double delta_dec)
Creates a resample_outgrid hdrl parameter object for a 2 dimensional interpolation,...
Definition hdrl_resample.c:2757
hdrl_resample_parameter_create_outgrid2D_userdef
hdrl_parameter * hdrl_resample_parameter_create_outgrid2D_userdef(const double delta_ra, const double delta_dec, const double ra_min, const double ra_max, const double dec_min, const double dec_max, const double fieldmargin)
Creates a resample_outgrid hdrl parameter object for a 2 dimensional interpolation,...
Definition hdrl_resample.c:2879
hdrl_resample_parameter_outgrid_verify
cpl_error_code hdrl_resample_parameter_outgrid_verify(const hdrl_parameter *hp)
verify parameters have proper values
Definition hdrl_resample.c:3292
hdrl_resample_compute
hdrl_resample_result * hdrl_resample_compute(const cpl_table *ResTable, hdrl_parameter *method, hdrl_parameter *outputgrid, const cpl_wcs *wcs)
High level resampling function.
Definition hdrl_resample.c:735
hdrl_resample_result_delete
void hdrl_resample_result_delete(hdrl_resample_result *aCube)
Deallocates the memory associated to a hdrl_resample_result object.
Definition hdrl_resample.c:848
hdrl_resample_image_to_table
cpl_table * hdrl_resample_image_to_table(const hdrl_image *hima, const cpl_wcs *wcs)
Convert a hdrl image into a cpl table that can be given as input to hdrl_resample_compute()
Definition hdrl_resample.c:2605
hdrl_resample_imagelist_to_table
cpl_table * hdrl_resample_imagelist_to_table(const hdrl_imagelist *himlist, const cpl_wcs *wcs)
Convert a hdrl imagelist into a cpl table that can be given as input to hdrl_resample_compute().
Definition hdrl_resample.c:2635
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2