CR2RE Pipeline Reference Manual 1.7.0
cr2res_calib-test.c
1/*
2 * This file is part of the CR2RES Pipeline
3 * Copyright (C) 2002,2003 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 02111-1307 USA
18 */
19
20
21#ifdef HAVE_CONFIG_H
22#include <config.h>
23#endif
24
25/*-----------------------------------------------------------------------------
26 Includes
27 -----------------------------------------------------------------------------*/
28
29#include <stdlib.h>
30#include <string.h>
31#include <math.h>
32#include <cpl.h>
33#include <hdrl.h>
34#include <cr2res_dfs.h>
35#include <cr2res_calib.h>
36#include <cr2res_io.h>
37#include "cr2res_pfits.h"
38#include "cr2res_utils.h"
39
40#define MODE_FLAT 0
41#define MODE_DARK 1
42#define MODE_BPM 2
43#define MODE_DETLIN 3
44
45#define pow2(x) (x) * (x)
46#define detlin(x, a, b, c) (x) * ((a) + (x) * ((b) + (x) * (c)))
47#define deterr(x, a, b, c, sx, sa, sb, sc) sqrt(detlin(pow2(x), pow2(sa), pow2(sb), pow2(sc)) + pow2(sx) * pow2((a) + 2 * (b) * (x) + 3 * (c) * pow2(x)))
48
49#define localdir "." //getenv("srcdir")
50
51static hdrl_image * cr2res_create_hdrl(int nx, int ny, double value, double error);
52
53static void test_cr2res_calib_image(void);
54static void test_cr2res_calib_cosmic(void);
55static void test_cr2res_calib_flat(void);
56static void test_cr2res_calib_dark(void);
57static void test_cr2res_calib_bpm(void);
58static void test_cr2res_calib_detlin(void);
59
60
61static void create_empty_fits()
62{
63 // Many tests require a empty fits file, since crires code checks if a file exists
64 // when adding a filename to a frame
65 char *my_path = cpl_sprintf("%s/TEST_empty.fits", localdir);
66 cpl_propertylist * list = cpl_propertylist_new();
67 cpl_propertylist_append_string(list, CR2RES_HEADER_DRS_TYPE, "DEBUG");
68 cpl_propertylist_save(list, my_path, CPL_IO_CREATE);
69 cpl_propertylist_delete(list);
70 cpl_free(my_path);
71}
72
73/*----------------------------------------------------------------------------*/
77/*----------------------------------------------------------------------------*/
78static void save_hdrl(char * filename, hdrl_image * hdrl, int mode, double dit)
79{
80 // Empty structures needed for master flat frame, but not for the test
81 // Need an empty fits file with header for the framesets
82 cpl_frame * empty = cpl_frame_new();
83
84 char *my_path = cpl_sprintf("%s/TEST_empty.fits", localdir);
85 cpl_frame_set_filename(empty, my_path);
86 cpl_frame_set_tag(empty, "DEBUG");
87 cpl_frame_set_group(empty, CPL_FRAME_GROUP_CALIB);
88
89 cpl_parameterlist * parlist = cpl_parameterlist_new();
90 cpl_propertylist * ext1 = cpl_propertylist_new();
91 cpl_propertylist_append_double(ext1, CR2RES_HEADER_DIT, dit);
92 cpl_propertylist * ext[] = {ext1, NULL, NULL};
93
94
95 cpl_frameset * all = cpl_frameset_new();
96 cpl_frameset * in = cpl_frameset_new();
97 cpl_frameset_insert(all, empty);
98 empty = cpl_frame_duplicate(empty);
99 cpl_frameset_insert(in, empty);
100
101 hdrl_image * list[] = {hdrl, NULL, NULL};
102
103 if (mode == MODE_FLAT)
104 cr2res_io_save_MASTER_FLAT(filename, all, in, parlist, list, ext1, ext, CR2RES_CAL_FLAT_MASTER_PROCATG, "debug");
105 if (mode == MODE_DARK)
106 cr2res_io_save_MASTER_DARK(filename, all, in, parlist, list, ext1, ext, CR2RES_CAL_DARK_MASTER_PROCATG, "debug");
107 if (mode == MODE_BPM)
108 {
109 cpl_image * list2[] = {hdrl_image_get_image(hdrl), NULL, NULL};
110 cr2res_io_save_BPM(filename, all, in, parlist, list2, ext1, ext, CR2RES_CAL_FLAT_BPM_PROCATG, "debug");
111 }
112 if (mode == MODE_DETLIN){
113 hdrl_imagelist * list3 = hdrl_imagelist_new();
114 hdrl_imagelist_set(list3, hdrl, 0);
115 hdrl_imagelist_set(list3, hdrl, 1);
116 hdrl_imagelist_set(list3, hdrl, 2);
117
118 hdrl_imagelist * list4[] = {list3, NULL, NULL};
119
120 cr2res_io_save_DETLIN_COEFFS(filename, all, in, parlist, list4, ext1, ext, CR2RES_CAL_DETLIN_COEFFS_PROCATG, "debug");
121
122 // hdrl_imagelist_unset(list3, 0);
123 // hdrl_imagelist_unset(list3, 1);
124 // hdrl_imagelist_unset(list3, 2);
125 hdrl_imagelist_delete(list3);
126 }
127 cpl_frameset_delete(all);
128 cpl_frameset_delete(in);
129 cpl_free(my_path) ;
130 cpl_parameterlist_delete(parlist);
131 cpl_propertylist_delete(ext1);
132}
133
134/*----------------------------------------------------------------------------*/
138/*----------------------------------------------------------------------------*/
139static cpl_frame * create_master_flat(char * filename, int nx, int ny, double value, double error, cpl_mask ** bpm)
140{
141 cpl_frame * out = cpl_frame_new();
142 cpl_frame_set_filename(out, filename);
143 cpl_frame_set_tag(out, "FLAT");
144 cpl_frame_set_group(out, CPL_FRAME_GROUP_CALIB);
145
146 hdrl_image * hdrl = cr2res_create_hdrl(nx, ny, value, error);
147
148 if (bpm != NULL){
149 hdrl_image_reject(hdrl, 1, 1);
150 *bpm = cpl_mask_duplicate(hdrl_image_get_mask(hdrl));
151 }
152
153 save_hdrl(filename, hdrl, MODE_FLAT, 0);
154 hdrl_image_delete(hdrl);
155
156 return out;
157}
158
159/*----------------------------------------------------------------------------*/
163/*----------------------------------------------------------------------------*/
164static cpl_frame * create_master_dark(char * filename, int nx, int ny, double value, double error, double dit, cpl_mask ** bpm)
165{
166
167 cpl_frame * out = cpl_frame_new();
168 cpl_frame_set_filename(out, filename);
169 cpl_frame_set_tag(out, "DARK");
170 cpl_frame_set_group(out, CPL_FRAME_GROUP_CALIB);
171
172 hdrl_image * hdrl = cr2res_create_hdrl(nx, ny, value, error);
173
174 if (bpm != NULL){
175 hdrl_image_reject(hdrl, 1, 1);
176 *bpm = cpl_mask_duplicate(hdrl_image_get_mask(hdrl));
177 }
178
179 save_hdrl(filename, hdrl, MODE_DARK, dit);
180 hdrl_image_delete(hdrl);
181
182 return out;
183}
184
185/*----------------------------------------------------------------------------*/
189/*----------------------------------------------------------------------------*/
190static cpl_frame * create_bpm(char * filename, int nx, int ny, double value)
191{
192 cpl_frame * out = cpl_frame_new();
193 cpl_frame_set_filename(out, filename);
194 cpl_frame_set_tag(out, "BPM");
195 cpl_frame_set_group(out, CPL_FRAME_GROUP_CALIB);
196
197 hdrl_image * hdrl = cr2res_create_hdrl(nx, ny, value, 0);
198 save_hdrl(filename, hdrl, MODE_BPM, 0);
199 hdrl_image_delete(hdrl);
200
201 return out;
202}
203
204/*----------------------------------------------------------------------------*/
208/*----------------------------------------------------------------------------*/
209static cpl_frame * create_detlin(char * filename, hdrl_image * a, hdrl_image * b, hdrl_image * c)
210{
211 char *my_path = cpl_sprintf("%s/TEST_empty.fits", localdir);
212
213 // Empty structures needed for master flat frame, but not for the test
214 // Need an empty fits file with header for the framesets
215 cpl_frame * empty = cpl_frame_new();
216 cpl_frame_set_filename(empty, my_path);
217 cpl_frame_set_tag(empty, "DEBUG");
218 cpl_frame_set_group(empty, CPL_FRAME_GROUP_CALIB);
219
220 cpl_parameterlist * parlist = cpl_parameterlist_new();
221 cpl_propertylist * ext1 = cpl_propertylist_new();
222 cpl_propertylist * ext[] = {ext1, NULL, NULL};
223
224 cpl_frameset * all = cpl_frameset_new();
225 cpl_frameset * in = cpl_frameset_new();
226 cpl_frameset_insert(all, empty);
227 empty = cpl_frame_duplicate(empty);
228 cpl_frameset_insert(in, empty);
229
230 hdrl_imagelist * list3 = hdrl_imagelist_new();
231 hdrl_imagelist_set(list3, a, 0);
232 hdrl_imagelist_set(list3, b, 1);
233 hdrl_imagelist_set(list3, c, 2);
234
235 hdrl_imagelist * list4[] = {list3, NULL, NULL};
236
237 cr2res_io_save_DETLIN_COEFFS(filename, all, in, parlist, list4, ext1, ext, CR2RES_CAL_DETLIN_COEFFS_PROCATG, "debug");
238
239 hdrl_imagelist_unset(list3, 2);
240 hdrl_imagelist_unset(list3, 1);
241 hdrl_imagelist_unset(list3, 0);
242
243 hdrl_imagelist_delete(list3);
244
245 cpl_frameset_delete(all);
246 cpl_frameset_delete(in);
247 cpl_parameterlist_delete(parlist);
248 cpl_propertylist_delete(ext1);
249
250 cpl_frame * out = cpl_frame_new();
251 cpl_frame_set_filename(out, filename);
252 cpl_frame_set_tag(out, "DETLIN");
253 cpl_frame_set_group(out, CPL_FRAME_GROUP_CALIB);
254
255 cpl_free(my_path);
256
257 return out;
258}
259
260/*----------------------------------------------------------------------------*/
264/*----------------------------------------------------------------------------*/
265static hdrl_image * cr2res_create_hdrl(int nx, int ny, double value, double error)
266{
267 hdrl_image * out = hdrl_image_new(nx, ny);
268 hdrl_value hv;
269 hv.data = value;
270 hv.error = error;
271 hdrl_image_add_scalar(out, hv);
272 return out;
273}
274
275static double calc_read_noise(double dit, int ndit, int chip)
276{
277 double min_dit = 1.427;
278 double lim_dit = 50.0;
279 double min_rn;
280 if (chip ==1){
281 min_rn = 11.0;
282 } else {
283 min_rn = 12.0;
284 }
285 double lim_rn = 6.0;
286 if(dit <= min_dit){
287 return min_rn/CR2RES_GAIN_CHIP1;
288 } else if(dit >= lim_dit) {
289 return lim_rn/CR2RES_GAIN_CHIP1;
290 } else {
291 double rn = min_rn + (lim_rn-min_rn)*(dit-min_dit)/(lim_dit-min_dit);
292 return rn/CR2RES_GAIN_CHIP1;
293 }
294}
295
296/* Expected cr2res_add_shotnoise error for chip 1 */
297static double expected_noise(double err, double value, double dit, int ndit,
298 int ndsamples)
299{
300 double n = ndsamples ;
301 double alpha = ndsamples < 2 ? 1.0 :
302 6.0 * (n * n + 1.0) / (5.0 * n * (n + 1.0)) ;
303 double rn = calc_read_noise(dit, 1, 1) ;
304 return sqrt(err * err + alpha * (value > 0.0 ? value : 0.0) /
305 (CR2RES_GAIN_CHIP1 * ndit) + rn * rn / ndit) ;
306}
307
308/* Applies expected_noise() with NDIT=1, NDSAMPLES=2 to every pixel of im */
309static void add_expected_noise(hdrl_image * im, double dit)
310{
311 cpl_image * err = hdrl_image_get_error(im) ;
312 const cpl_image * val = hdrl_image_get_image_const(im) ;
313 cpl_size i, j ;
314 int badpix ;
315
316 for (j = 1; j <= cpl_image_get_size_y(err); j++) {
317 for (i = 1; i <= cpl_image_get_size_x(err); i++) {
318 cpl_image_set(err, i, j, expected_noise(
319 cpl_image_get(err, i, j, &badpix),
320 cpl_image_get(val, i, j, &badpix), dit, 1, 2)) ;
321 }
322 }
323}
324
325static void test_cr2res_calib_image()
326{
327 int nx = 5;
328 int ny = 5;
329 //int badpix;
330
331 double img_value = 100;
332 double img_error = 1;
333
334 hdrl_image * in = cr2res_create_hdrl(nx, ny, img_value, img_error);
335 int chip = 1;
336
337 char *my_path1 = cpl_sprintf("%s/TEST_master_flat.fits", localdir);
338 char *my_path2 = cpl_sprintf("%s/TEST_master_dark.fits", localdir);
339 char *my_path3 = cpl_sprintf("%s/TEST_bpm.fits", localdir);
340 cpl_frame * flat = create_master_flat(my_path1, nx, ny, 1, 0, NULL);
341 cpl_frame * dark = create_master_dark(my_path2, nx, ny, 10, 1, 10, NULL);
342 cpl_frame * bpm = create_bpm(my_path3, nx, ny, 0);
343 cpl_frame * detlin = NULL;
344 double dit = 10;
345
346 hdrl_image * out;
347 hdrl_image * cmp;
348
349 // NULL input / output
350 out = cr2res_calib_image(NULL, chip, 0, 0, 1, 0, NULL, NULL, NULL, NULL,dit,1,2);
351 cpl_test_null(out);
352
353 out = cr2res_calib_image(in, 0, 0, 0, 1, 0, NULL, NULL, NULL, NULL, dit, 1, 2);
354 cpl_test_null(out);
355
356 out = cr2res_calib_image(in, CR2RES_NB_DETECTORS + 1, 0, 0, 1, 0,
357 NULL, NULL, NULL, NULL, dit, 1, 2);
358 cpl_test_null(out);
359
360 // No correction
361 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, NULL, NULL, NULL, NULL, dit, 1, 2);
362 cpl_test_nonnull(out);
363 cpl_test_image_abs(hdrl_image_get_image(in), hdrl_image_get_image(out), DBL_EPSILON);
364
365 // shot and read noise are added in quadrature
366 cmp = hdrl_image_duplicate(in);
367 add_expected_noise(cmp, dit);
368
369 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
370 hdrl_image_delete(out);
371 hdrl_image_delete(cmp);
372
373 // NDIT and up-the-ramp NDSAMPLES
374 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, NULL, NULL, NULL, NULL, dit, 4, 33);
375 cpl_test_nonnull(out);
376 cmp = cr2res_create_hdrl(nx, ny, img_value,
377 expected_noise(img_error, img_value, dit, 4, 33));
378 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
379 hdrl_image_delete(out);
380 hdrl_image_delete(cmp);
381
382 // Negative signal adds no shot noise
383 hdrl_image * neg = cr2res_create_hdrl(nx, ny, -img_value, img_error);
384 out = cr2res_calib_image(neg, chip, 0, 0, 1, 0, NULL, NULL, NULL, NULL, dit, 1, 2);
385 cpl_test_nonnull(out);
386 cmp = cr2res_create_hdrl(nx, ny, -img_value,
387 expected_noise(img_error, -img_value, dit, 1, 2));
388 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
389 hdrl_image_delete(out);
390 hdrl_image_delete(cmp);
391 hdrl_image_delete(neg);
392
393
394 hdrl_image_delete(in);
395 cpl_frame_delete(flat);
396 cpl_frame_delete(dark);
397 cpl_frame_delete(bpm);
398 cpl_frame_delete(detlin);
399 cpl_free(my_path1) ;
400 cpl_free(my_path2) ;
401 cpl_free(my_path3) ;
402}
403
404/*----------------------------------------------------------------------------*/
408/*----------------------------------------------------------------------------*/
409static void test_cr2res_calib_cosmic()
410{
411 int nx = 5;
412 int ny = 5;
413 //int badpix;
414
415 double img_value = 100;
416 double img_error = 1;
417 //double out_value, out_error;
418 //double flat_value, flat_error;
419
420 hdrl_image * in, * out;// * cmp;
421 //cpl_frame * flat;
422 int chip = 1;
423 //int cosmics_corr;
424 double dit = 10;
425
426 in = cr2res_create_hdrl(nx, ny, img_value, img_error);
427
428 // Case 1: No Cosmic Correction, i.e. nothing happens
429 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, NULL, NULL, NULL, NULL, dit, 1, 2);
430 cpl_test_image_abs(hdrl_image_get_image(in), hdrl_image_get_image(out), DBL_EPSILON);
431 add_expected_noise(in, dit);
432
433 cpl_test_image_abs(hdrl_image_get_error(in), hdrl_image_get_error(out), 1e-12);
434
435 // Case 2: Cosmic Correction
436 // TODO when cosmic correction is implemented
437
438 hdrl_image_delete(in);
439 hdrl_image_delete(out);
440
441}
442
443
444/*----------------------------------------------------------------------------*/
448/*----------------------------------------------------------------------------*/
449static void test_cr2res_calib_flat()
450{
451 int nx = 5;
452 int ny = 5;
453 //int badpix;
454
455 double img_value = 100;
456 double img_error = 1;
457 double tmp_error = 0;
458 double out_value, out_error;
459 double flat_value, flat_error;
460
461 hdrl_image * in, * out, * cmp;
462 //cpl_mask * bpm;
463 cpl_frame * flat;
464 int chip = 1;
465 double dit = 10;
466 //cpl_error_code error;
467
468 in = cr2res_create_hdrl(nx, ny, img_value, img_error);
469
470 // Case 1: Flat is just 1 and no error, i.e. no change
471 flat_value = 1;
472 flat_error = 0;
473 char *my_path1 = cpl_sprintf("%s/TEST_master_flat.fits", localdir);
474 flat = create_master_flat(my_path1, nx, ny, flat_value, flat_error, NULL);
475 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, flat, NULL, NULL, NULL, dit, 1, 2);
476
477 out_value = img_value / flat_value;
478 out_error = sqrt( pow(img_error / flat_value, 2) + pow(img_value * flat_error/ (flat_value * flat_value), 2));
479 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
480 add_expected_noise(cmp, dit);
481
482 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
483 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
484
485 hdrl_image_delete(out);
486 hdrl_image_delete(cmp);
487 cpl_frame_delete(flat);
488
489 // Case 2: Flat is constant == 2, error 1
490 flat_value = 2;
491 flat_error = 1;
492 flat = create_master_flat(my_path1, nx, ny, flat_value, flat_error, NULL);
493 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, flat, NULL, NULL, NULL, dit, 1, 2);
494
495 out_value = img_value / flat_value;
496 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
497
498 out_error = sqrt( pow(tmp_error / flat_value, 2) + pow(img_value * flat_error/ (flat_value * flat_value), 2));
499 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
500
501 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
502 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
503
504 hdrl_image_delete(out);
505 hdrl_image_delete(cmp);
506 cpl_frame_delete(flat);
507
508 // Case 3: Flat is 0, i.e. all pixels become bad
509 flat_value = 0;
510 flat_error = 1;
511 flat = create_master_flat(my_path1, nx, ny, flat_value, flat_error, NULL);
512 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, flat, NULL, NULL, NULL, dit, 1, 2);
513
514 cpl_test_eq(nx * ny, hdrl_image_count_rejected(out));
515
516 hdrl_image_delete(out);
517 cpl_frame_delete(flat);
518
519 // Case 4: Flat file does not exist
520 flat = cpl_frame_new();
521 cpl_frame_set_filename(flat, "TEST_tobeornottobe.fits");
522 cpl_frame_set_tag(flat, "FLAT");
523 cpl_frame_set_group(flat, CPL_FRAME_GROUP_CALIB);
524
525 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, flat, NULL, NULL, NULL, dit, 1, 2);
526 cpl_test_error(CPL_ERROR_FILE_IO);
527 hdrl_image_delete(out);
528 cpl_frame_delete(flat);
529
530 // Case 5: image is in a wrong group
531 flat = create_master_dark(my_path1, nx, ny, 1, 0, 10, NULL);
532 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, flat, NULL, NULL, NULL, dit, 1, 2);
533 cpl_test_null(out);
534 cpl_frame_delete(flat);
535
536 // Case 6: No Filename set
537 flat = cpl_frame_new();
538 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, flat, NULL, NULL, NULL, dit, 1, 2);
539 cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
540 cpl_test_null(out);
541 cpl_frame_delete(flat);
542
543 // Case 7: DataFile is empty, i.e. only header
544 flat = cpl_frame_new();
545 char *my_path2 = cpl_sprintf("%s/TEST_empty.fits", localdir);
546 cpl_frame_set_filename(flat, my_path2);
547 out = cr2res_calib_image(in, chip, 0, 0, 1, 0, flat, NULL, NULL, NULL, dit, 1, 2);
548 cpl_test_null(out);
549 cpl_frame_delete(flat);
550
551 // // Case 8: BPM is set
552 // flat_value = 1;
553 // flat_error = 1;
554 // flat = create_master_flat("TEST_master_flat.fits", nx, ny, flat_value, flat_error, &bpm);
555 // out = cr2res_calib_image(in, chip, 0, 0, 0, flat, NULL, NULL, NULL, dit, 1);
556
557 // cpl_test_eq_mask(bpm, hdrl_image_get_mask(out));
558
559 // hdrl_image_delete(out);
560 // cpl_frame_delete(flat);
561 // cpl_mask_delete(bpm);
562
563
564 cpl_free(my_path1) ;
565 cpl_free(my_path2) ;
566 hdrl_image_delete(in);
567}
568
569static void test_cr2res_calib_dark()
570{
571 int nx = 5;
572 int ny = 5;
573 //int badpix;
574
575 double img_value = 100;
576 double img_error = 1;
577 double tmp_error = 0;
578 double out_value, out_error;
579 double dark_value, dark_error, dark_dit;
580
581 hdrl_image * in, * out, * cmp;
582 //cpl_mask *bpm;
583 cpl_frame * dark;
584 int chip = 1;
585 double dit = 10;
586
587 in = cr2res_create_hdrl(nx, ny, img_value, img_error);
588 char *my_path = cpl_sprintf("%s/TEST_master_dark.fits", localdir);
589 char *my_path2 = cpl_sprintf("%s/TEST_empty.fits", localdir);
590
591
592 // Case 1: Dark is 0, dit is the same as the image, i.e. no correction
593 dark_value = 0;
594 dark_error = 0;
595 dark_dit = dit;
596 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
597 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
598
599 out_value = img_value - dit/dark_dit * dark_value;
600 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
601 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
602 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
603
604 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
605 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
606
607 hdrl_image_delete(out);
608 hdrl_image_delete(cmp);
609 cpl_frame_delete(dark);
610
611 // Case 2: Dark is 0, dit is different, i.e. still nothing
612 dark_value = 0;
613 dark_error = 0;
614 dark_dit = dit * 2.3;
615 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
616 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
617
618 out_value = img_value - dit/dark_dit * dark_value;
619 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
620 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
621 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
622
623 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
624 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
625
626 hdrl_image_delete(out);
627 hdrl_image_delete(cmp);
628 cpl_frame_delete(dark);
629
630 // Case 3: Dark is 0, dit is same, but dark has error now
631 dark_value = 0;
632 dark_error = 2;
633 dark_dit = dit;
634 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
635 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
636
637 out_value = img_value - dit/dark_dit * dark_value;
638 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
639 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
640 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
641
642 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
643 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
644
645 hdrl_image_delete(out);
646 hdrl_image_delete(cmp);
647 cpl_frame_delete(dark);
648
649 // Case 4: Dark is constant 10, dit is same, and has error
650 dark_value = 10;
651 dark_error = 2;
652 dark_dit = dit;
653 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
654 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
655
656 out_value = img_value - dit/dark_dit * dark_value;
657 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
658 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
659 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
660
661 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
662 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
663
664 hdrl_image_delete(out);
665 hdrl_image_delete(cmp);
666 cpl_frame_delete(dark);
667
668 // Case 5: Dark is constant 10, dit is different, and has error
669 dark_value = 10;
670 dark_error = 2;
671 dark_dit = dit * 2.3;
672 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
673 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
674
675 out_value = img_value - dit/dark_dit * dark_value;
676 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
677 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
678 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
679
680 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
681 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
682
683 hdrl_image_delete(out);
684 hdrl_image_delete(cmp);
685 cpl_frame_delete(dark);
686
687 // Case 6: Dark is negative, dit positive
688 dark_value = -10;
689 dark_error = 2;
690 dark_dit = dit * 2.3;
691 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
692 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
693
694 out_value = img_value - dit/dark_dit * dark_value;
695 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
696 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
697 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
698
699 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
700 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
701
702 hdrl_image_delete(out);
703 hdrl_image_delete(cmp);
704 cpl_frame_delete(dark);
705
706 // Case 7: Dark is positive, dit negative
707 dark_value = 10;
708 dark_error = 2;
709 dark_dit = -dit * 2.3;
710 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
711 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
712
713 out_value = img_value - dit/dark_dit * dark_value;
714 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
715 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
716 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
717
718 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
719 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
720
721 hdrl_image_delete(out);
722 hdrl_image_delete(cmp);
723 cpl_frame_delete(dark);
724
725 // Case 8: Dark and dit negative
726 dark_value = -10;
727 dark_error = 2;
728 dark_dit = -dit * 2.3;
729 dark = create_master_dark(my_path, nx, ny, dark_value, dark_error, dark_dit, NULL);
730 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
731
732 out_value = img_value - dit/dark_dit * dark_value;
733 tmp_error = expected_noise(img_error, img_value, dit, 1, 2);
734 out_error = sqrt(pow(tmp_error, 2) + pow(dit/dark_dit * dark_error, 2));
735 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
736
737 cpl_test_image_abs(hdrl_image_get_image(cmp), hdrl_image_get_image(out), DBL_EPSILON);
738 cpl_test_image_abs(hdrl_image_get_error(cmp), hdrl_image_get_error(out), 1e-12);
739
740 hdrl_image_delete(out);
741 hdrl_image_delete(cmp);
742 cpl_frame_delete(dark);
743
744 // Case 9: Dark file does not exist
745 dark = cpl_frame_new();
746 cpl_frame_set_filename(dark, "TEST_tobeornottobe.fits");
747 cpl_frame_set_tag(dark, "DARK");
748 cpl_frame_set_group(dark, CPL_FRAME_GROUP_CALIB);
749
750 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
751 cpl_test_error(CPL_ERROR_FILE_IO);
752 cpl_test_null(out);
753 cpl_frame_delete(dark);
754
755 // Case 10: image is in a wrong group
756 dark = create_master_flat(my_path, nx, ny, 1, 0, NULL);
757 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
758 cpl_test_null(out);
759 cpl_frame_delete(dark);
760
761 // Case 11: No Filename set
762 dark = cpl_frame_new();
763 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
764 cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
765 cpl_test_null(out);
766 cpl_frame_delete(dark);
767
768
769
770 // Case 12: DataFile is empty, i.e. only header
771 dark = cpl_frame_new();
772 cpl_frame_set_filename(dark, my_path2);
773 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1, 2);
774 cpl_test_null(out);
775 cpl_frame_delete(dark);
776
777 // // Case 13: Check that bpm of dark is passed on to out
778 // dark_value = 10;
779 // dark_error = 0;
780 // dark_dit = 10;
781 // dark = create_master_dark("TEST_master_dark.fits", nx, ny, dark_value, dark_error, dark_dit, &bpm);
782 // out = cr2res_calib_image(in, chip, 0, 0, 0, NULL, dark, NULL, NULL, dit, 1);
783
784 // cpl_test_eq_mask(bpm, hdrl_image_get_mask(out));
785
786 // cpl_mask_delete(bpm);
787 // hdrl_image_delete(out);
788 // cpl_frame_delete(dark);
789
790 cpl_free(my_path);
791 cpl_free(my_path2);
792
793 hdrl_image_delete(in);
794}
795
796static void test_cr2res_calib_bpm()
797{
798 int nx = 5;
799 int ny = 5;
800 int badpix;
801
802 double img_value = 100;
803 double img_error = 1;
804 //double out_value, out_error;
805 //double bpm_value;
806 hdrl_image * hdrl;
807 hdrl_value value;
808
809 hdrl_image * in, * out, * cmp;
810 cpl_frame * bpm;
811 int chip = 1;
812 double dit = 10;
813
814 in = cr2res_create_hdrl(nx, ny, img_value, img_error);
815 value.data = 10; value.error = img_error;
816 hdrl_image_set_pixel(in, 2, 2, value);
817
818 char *my_path = cpl_sprintf("%s/TEST_master_bpm.fits", localdir);
819 char *my_path2 = cpl_sprintf("%s/TEST_empty.fits", localdir);
820
821
822 // Case 1: BPM is all 0, i.e. all good pixels
823 bpm = create_bpm(my_path, nx, ny, 0);
824 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
825 cmp = hdrl_image_duplicate(in);
826 add_expected_noise(cmp, dit);
827
828 cpl_test_image_abs(hdrl_image_get_image(out), hdrl_image_get_image(cmp), DBL_EPSILON);
829 cpl_test_image_abs(hdrl_image_get_error(out), hdrl_image_get_error(cmp), 1e-12);
830 cpl_frame_delete(bpm);
831 hdrl_image_delete(out);
832 hdrl_image_delete(cmp);
833
834 // Case 2: BPM is all 1, i.e all pixels are bad
835 bpm = create_bpm(my_path, nx, ny, 1);
836 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
837 cpl_test_nonnull(out);
838 cpl_test_eq(hdrl_image_count_rejected(out), nx * ny);
839 cpl_frame_delete(bpm);
840 hdrl_image_delete(out);
841
842 // Case 3: One Bad Pixel
843 bpm = cpl_frame_new();
844 cpl_frame_set_filename(bpm, my_path);
845 cpl_frame_set_tag(bpm, "BPM");
846 cpl_frame_set_group(bpm, CPL_FRAME_GROUP_CALIB);
847
848 hdrl = cr2res_create_hdrl(nx, ny, 0, 0);
849 value.data = 1;
850 hdrl_image_set_pixel(hdrl, 2, 2, value);
851 save_hdrl(my_path, hdrl, MODE_BPM, 0);
852 hdrl_image_delete(hdrl);
853
854 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
855 cmp = hdrl_image_duplicate(in);
856 add_expected_noise(cmp, dit);
857
858 cpl_test_eq(hdrl_image_count_rejected(out), 1);
859 cpl_test_image_abs(hdrl_image_get_image(out), hdrl_image_get_image(cmp), DBL_EPSILON);
860 // cpl_test_image_abs(hdrl_image_get_error(out), hdrl_image_get_error(cmp), DBL_EPSILON);
861
862 for (size_t i = 1; i <= (size_t) nx; i++)
863 {
864 for (size_t j = 1; j <= (size_t) ny; j++)
865 {
866 if (!hdrl_image_is_rejected(out, i, j)){
867 cpl_test_abs(cpl_image_get(hdrl_image_get_error(out), i, j, &badpix),
868 cpl_image_get(hdrl_image_get_error(cmp), i, j, &badpix),
869 1e-12);
870 }
871 }
872
873 }
874
875 cpl_frame_delete(bpm);
876 hdrl_image_delete(out);
877 hdrl_image_delete(cmp);
878
879 // Case 4: Only one good pixel
880 bpm = cpl_frame_new();
881 cpl_frame_set_filename(bpm, my_path);
882 cpl_frame_set_tag(bpm, "BPM");
883 cpl_frame_set_group(bpm, CPL_FRAME_GROUP_CALIB);
884
885 hdrl = cr2res_create_hdrl(nx, ny, 1, 0);
886 value.data = 0;
887 hdrl_image_set_pixel(hdrl, 2, 2, value);
888 save_hdrl(my_path, hdrl, MODE_BPM, 0);
889 hdrl_image_delete(hdrl);
890
891 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
892 cmp = hdrl_image_duplicate(in);
893 add_expected_noise(cmp, dit);
894
895 cpl_test_eq(hdrl_image_count_rejected(out), nx*ny - 1);
896 cpl_test_image_abs(hdrl_image_get_image(out), hdrl_image_get_image(cmp), DBL_EPSILON);
897 // cpl_test_image_abs(hdrl_image_get_error(out), hdrl_image_get_error(cmp), DBL_EPSILON);
898
899 for (size_t i = 1; i <= (size_t) nx; i++)
900 {
901 for (size_t j = 1; j <= (size_t) ny; j++)
902 {
903 if (!hdrl_image_is_rejected(out, i, j)){
904 cpl_test_abs(cpl_image_get(hdrl_image_get_error(out), i, j, &badpix),
905 cpl_image_get(hdrl_image_get_error(cmp), i, j, &badpix),
906 1e-12);
907 }
908 }
909
910 }
911
912
913 cpl_frame_delete(bpm);
914 hdrl_image_delete(out);
915 hdrl_image_delete(cmp);
916
917 // Case 5: BPM file does not exist
918 bpm = cpl_frame_new();
919 cpl_frame_set_filename(bpm, "TEST_tobeornottobe.fits");
920 cpl_frame_set_tag(bpm, "BPM");
921 cpl_frame_set_group(bpm, CPL_FRAME_GROUP_CALIB);
922
923 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
924 cpl_test_error(CPL_ERROR_FILE_IO);
925 cpl_test_null(out);
926 cpl_frame_delete(bpm);
927
928 // Case 6: image is in a wrong group
929 bpm = create_master_flat(my_path, nx, ny, 1, 0, NULL);
930 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
931 cpl_test_null(out);
932 cpl_frame_delete(bpm);
933
934 // Case 7: No Filename set
935 bpm = cpl_frame_new();
936 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
937 cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
938 cpl_test_null(out);
939 cpl_frame_delete(bpm);
940
941
942 // Case 8: DataFile is empty, i.e. only header
943 bpm = cpl_frame_new();
944 cpl_frame_set_filename(bpm, my_path2);
945 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, bpm, NULL, dit, 1, 2);
946 cpl_test_null(out);
947 cpl_frame_delete(bpm);
948
949 cpl_free(my_path);
950 cpl_free(my_path2);
951
952 hdrl_image_delete(in);
953}
954
955static void test_cr2res_calib_detlin()
956{
957 int nx = 5;
958 int ny = 5;
959 //int badpix;
960
961 double img_value = 20000;
962 double img_error = 0;
963 double out_value, out_error;
964 //double bpm_value;
965 //hdrl_image * hdrl;
966 //hdrl_value value;
967
968 hdrl_image * in, * out, * cmp;
969 hdrl_image * ima, * imb, * imc;
970 double a, b, c, sa, sb, sc;
971 double i, si;
972 cpl_frame * detlin;
973 int chip = 1;
974 double dit = 10;
975
976 in = cr2res_create_hdrl(nx, ny, img_value, img_error);
977 i = img_value; si = img_error;
978
979 char *my_path = cpl_sprintf("%s/TEST_master_detlin.fits", localdir);
980 char *my_path2 = cpl_sprintf("%s/TEST_empty.fits", localdir);
981
982 // DetLin images aren't actually images, but polynomial coefficients for each pixel
983 // The new image is the old multiplied by the
984 // evaluation of that polynomial
985 // new image = old_image * detlin_poly(old_image)
986
987 // Case 1: DetLin b = c = 0, a = 1, doesn't change anything
988 a = 1; c = 0; b = 0;
989 sa = sb = sc = 0;
990 ima = cr2res_create_hdrl(nx, ny, a, sa);
991 imb = cr2res_create_hdrl(nx, ny, b, sb);
992 imc = cr2res_create_hdrl(nx, ny, c, sc);
993
994 detlin = create_detlin(my_path, ima, imb, imc);
995 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
996 cmp = hdrl_image_duplicate(in);
997 add_expected_noise(cmp, dit);
998
999 cpl_test_image_abs(hdrl_image_get_image(out), hdrl_image_get_image(cmp), DBL_EPSILON);
1000 cpl_test_image_rel(hdrl_image_get_error(out), hdrl_image_get_error(cmp), 1e-12);
1001 hdrl_image_delete(out);
1002 cpl_frame_delete(detlin);
1003 hdrl_image_delete(ima);
1004 hdrl_image_delete(imb);
1005 hdrl_image_delete(imc);
1006 hdrl_image_delete(cmp);
1007
1008 // Case 2: DetLin a = 1, c = 0, b = 1, x + x**2
1009 a = b = 1; c = 0;
1010 sa = sb = sc = 0;
1011 ima = cr2res_create_hdrl(nx, ny, a, sa);
1012 imb = cr2res_create_hdrl(nx, ny, b, sb);
1013 imc = cr2res_create_hdrl(nx, ny, c, sc);
1014 detlin = create_detlin(my_path, ima, imb, imc);
1015 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1016 out_value = detlin(img_value, a, b, c);
1017 out_error = deterr(img_value, a, b, c, img_error, sa, sb, sc);
1018 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1019 add_expected_noise(cmp, dit);
1020
1021 cpl_test_image_abs(hdrl_image_get_image(out),
1022 hdrl_image_get_image(cmp), DBL_EPSILON);
1023 cpl_test_image_rel(hdrl_image_get_error(out),
1024 hdrl_image_get_error(cmp), 1e-12);
1025 hdrl_image_delete(out);
1026 hdrl_image_delete(cmp);
1027 cpl_frame_delete(detlin);
1028 hdrl_image_delete(ima);
1029 hdrl_image_delete(imb);
1030 hdrl_image_delete(imc);
1031
1032
1033 // Case 3: DetLin a = b = c = 0, everything is 0
1034 a = c = b = 0;
1035 sa = sb = sc = 0;
1036 ima = cr2res_create_hdrl(nx, ny, a, sa);
1037 imb = cr2res_create_hdrl(nx, ny, b, sb);
1038 imc = cr2res_create_hdrl(nx, ny, c, sc);
1039 detlin = create_detlin(my_path, ima, imb, imc);
1040 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1041 cmp = cr2res_create_hdrl(nx, ny, 0, 0);
1042 add_expected_noise(cmp, dit);
1043 // Shot noise is also 0
1044 cpl_test_image_abs(hdrl_image_get_image(out),
1045 hdrl_image_get_image(cmp), DBL_EPSILON);
1046 cpl_test_image_abs(hdrl_image_get_error(out),
1047 hdrl_image_get_error(cmp), 1e-12);
1048 hdrl_image_delete(out);
1049 hdrl_image_delete(cmp);
1050 cpl_frame_delete(detlin);
1051 hdrl_image_delete(ima);
1052 hdrl_image_delete(imb);
1053 hdrl_image_delete(imc);
1054
1055 // Case 4: DetLin a = 0.5, c = b = 0, divide all values by 2, including the error
1056 a = 0.5; b = c = 0;
1057 sa = sb = sc = 0;
1058 ima = cr2res_create_hdrl(nx, ny, a, sa);
1059 imb = cr2res_create_hdrl(nx, ny, b, sb);
1060 imc = cr2res_create_hdrl(nx, ny, c, sc);
1061 detlin = create_detlin(my_path, ima, imb, imc);
1062 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1063 out_value = img_value * a;
1064 out_error = img_error * a;
1065 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1066 add_expected_noise(cmp, dit);
1067
1068 cpl_test_image_abs(hdrl_image_get_image(out),
1069 hdrl_image_get_image(cmp), DBL_EPSILON);
1070 cpl_test_image_rel(hdrl_image_get_error(out),
1071 hdrl_image_get_error(cmp), 1e-12);
1072 hdrl_image_delete(out);
1073 hdrl_image_delete(cmp);
1074 cpl_frame_delete(detlin);
1075 hdrl_image_delete(ima);
1076 hdrl_image_delete(imb);
1077 hdrl_image_delete(imc);
1078
1079 // Case 5: DetLin a = 1, c = 0, b = 2
1080 // x = a + y * b
1081 a = 2; b = 0; c = 0;
1082 sa = sb = sc = 0;
1083 ima = cr2res_create_hdrl(nx, ny, a, sa);
1084 imb = cr2res_create_hdrl(nx, ny, b, sb);
1085 imc = cr2res_create_hdrl(nx, ny, c, sc);
1086 detlin = create_detlin(my_path, ima, imb, imc);
1087 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1088 out_value = detlin(i, a, b, c);
1089 out_error = deterr(i, a, b, c, si, sa, sb, sc);
1090 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1091 add_expected_noise(cmp, dit);
1092
1093 cpl_test_image_abs(hdrl_image_get_image(out),
1094 hdrl_image_get_image(cmp), DBL_EPSILON);
1095 cpl_test_image_rel(hdrl_image_get_error(out),
1096 hdrl_image_get_error(cmp), 1e-12);
1097 hdrl_image_delete(out);
1098 hdrl_image_delete(cmp);
1099 cpl_frame_delete(detlin);
1100 hdrl_image_delete(ima);
1101 hdrl_image_delete(imb);
1102 hdrl_image_delete(imc);
1103
1104 // Case 6: DetLin a = 1, c = 1, b = 2
1105 // x = a + y * b + y**2 * c
1106 a = c = 1; b = 2;
1107 sa = sb = sc = 0;
1108 ima = cr2res_create_hdrl(nx, ny, a, sa);
1109 imb = cr2res_create_hdrl(nx, ny, b, sb);
1110 imc = cr2res_create_hdrl(nx, ny, c, sc);
1111 detlin = create_detlin(my_path, ima, imb, imc);
1112 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1113 out_value = detlin(i, a, b, c);
1114 out_error = deterr(i, a, b, c, si, sa, sb, sc);
1115 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1116 add_expected_noise(cmp, dit);
1117
1118 cpl_test_image_abs(hdrl_image_get_image(out),
1119 hdrl_image_get_image(cmp), DBL_EPSILON);
1120 cpl_test_image_rel(hdrl_image_get_error(out),
1121 hdrl_image_get_error(cmp), 1e-12);
1122 hdrl_image_delete(out);
1123 hdrl_image_delete(cmp);
1124 cpl_frame_delete(detlin);
1125 hdrl_image_delete(ima);
1126 hdrl_image_delete(imb);
1127 hdrl_image_delete(imc);
1128
1129 // Case 7: DetLin a = 0, c = 0, b = 1, but now with error
1130 // x = b * y
1131 a = c = 0; b = 1;
1132 sa = sc = 0; sb = 1;
1133 ima = cr2res_create_hdrl(nx, ny, a, sa);
1134 imb = cr2res_create_hdrl(nx, ny, b, sb);
1135 imc = cr2res_create_hdrl(nx, ny, c, sc);
1136 detlin = create_detlin(my_path, ima, imb, imc);
1137 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1138 out_value = pow2(img_value);
1139 out_error = deterr(i, a, b, c, si, sa, sb, sc);
1140 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1141 add_expected_noise(cmp, dit);
1142
1143 cpl_test_image_abs(hdrl_image_get_image(out), hdrl_image_get_image(cmp), DBL_EPSILON);
1144 cpl_test_image_rel(hdrl_image_get_error(out), hdrl_image_get_error(cmp), 1e-12);
1145 hdrl_image_delete(out);
1146 hdrl_image_delete(cmp);
1147 cpl_frame_delete(detlin);
1148 hdrl_image_delete(ima);
1149 hdrl_image_delete(imb);
1150 hdrl_image_delete(imc);
1151
1152 // Case 8: DetLin a = 1, c = 0, b = 1, but now with error
1153 a = 1; c = 0; b = 1;
1154 sa = 1 ; sc = 0; sb = 1;
1155 ima = cr2res_create_hdrl(nx, ny, a, sa);
1156 imb = cr2res_create_hdrl(nx, ny, b, sb);
1157 imc = cr2res_create_hdrl(nx, ny, c, sc);
1158 detlin = create_detlin(my_path, ima, imb, imc);
1159 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1160 out_value = detlin(i, a, b, c);
1161 out_error = deterr(i, a, b, c, si, sa, sb, sc);
1162 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1163 add_expected_noise(cmp, dit);
1164
1165 cpl_test_image_abs(hdrl_image_get_image(out),
1166 hdrl_image_get_image(cmp), DBL_EPSILON);
1167 cpl_test_image_rel(hdrl_image_get_error(out),
1168 hdrl_image_get_error(cmp), 1e-12);
1169 hdrl_image_delete(out);
1170 hdrl_image_delete(cmp);
1171 cpl_frame_delete(detlin);
1172 hdrl_image_delete(ima);
1173 hdrl_image_delete(imb);
1174 hdrl_image_delete(imc);
1175
1176 // Case 9: DetLin a = 1, c = 1, b = 0, but now with error
1177 // x = a + c * y**2
1178 a = 1; c = 1; b = 0;
1179 sa = 1 ; sc = 1; sb = 0;
1180 i = img_value; si = img_error;
1181 ima = cr2res_create_hdrl(nx, ny, a, sa);
1182 imb = cr2res_create_hdrl(nx, ny, b, sb);
1183 imc = cr2res_create_hdrl(nx, ny, c, sc);
1184 detlin = create_detlin(my_path, ima, imb, imc);
1185 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1186 out_value = detlin(i, a, b, c);
1187 out_error = deterr(i, a, b, c, si, sa, sb, sc);
1188 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1189 add_expected_noise(cmp, dit);
1190
1191 cpl_test_image_abs(hdrl_image_get_image(out), hdrl_image_get_image(cmp), DBL_EPSILON);
1192 cpl_test_image_rel(hdrl_image_get_error(out), hdrl_image_get_error(cmp), 1e-12);
1193 hdrl_image_delete(out);
1194 hdrl_image_delete(cmp);
1195 cpl_frame_delete(detlin);
1196 hdrl_image_delete(ima);
1197 hdrl_image_delete(imb);
1198 hdrl_image_delete(imc);
1199
1200 // Case 10: DetLin a = 1, c = 1, b = 1, but now with error
1201 // x = a + b * y + c * y**2
1202 a = 1; c = 1; b = 1;
1203 sa = 1 ; sc = 1; sb = 1;
1204 i = img_value; si = img_error;
1205 ima = cr2res_create_hdrl(nx, ny, a, sa);
1206 imb = cr2res_create_hdrl(nx, ny, b, sb);
1207 imc = cr2res_create_hdrl(nx, ny, c, sc);
1208 detlin = create_detlin(my_path, ima, imb, imc);
1209 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1210 out_value = detlin(i, a, b, c);
1211 out_error = deterr(i, a, b, c, si, sa, sb, sc);
1212 cmp = cr2res_create_hdrl(nx, ny, out_value, out_error);
1213 add_expected_noise(cmp, dit);
1214
1215 cpl_test_image_abs(hdrl_image_get_image(out), hdrl_image_get_image(cmp), DBL_EPSILON);
1216 cpl_test_image_rel(hdrl_image_get_error(out), hdrl_image_get_error(cmp), 1e-12);
1217 hdrl_image_delete(out);
1218 hdrl_image_delete(cmp);
1219 cpl_frame_delete(detlin);
1220 hdrl_image_delete(ima);
1221 hdrl_image_delete(imb);
1222 hdrl_image_delete(imc);
1223
1224 // Case 11: No Filename set
1225 detlin = cpl_frame_new();
1226 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1227 cpl_test_error(CPL_ERROR_DATA_NOT_FOUND);
1228 cpl_test_null(out);
1229 cpl_frame_delete(detlin);
1230
1231
1232 // Case 12: DataFile is empty, i.e. only header
1233 detlin = cpl_frame_new();
1234 cpl_frame_set_filename(detlin, my_path2);
1235 out = cr2res_calib_image(in, chip, 0, 0, 0, 0, NULL, NULL, NULL, detlin, dit, 1, 2);
1236 cpl_test_null(out);
1237 cpl_frame_delete(detlin);
1238
1239 cpl_free(my_path);
1240 cpl_free(my_path2);
1241
1242 hdrl_image_delete(in);
1243
1244}
1245
1246/*----------------------------------------------------------------------------*/
1250/*----------------------------------------------------------------------------*/
1251int main(void)
1252{
1253 cpl_error_reset();
1254 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_DEBUG);
1255
1256 create_empty_fits();
1257
1258 test_cr2res_calib_image();
1259 test_cr2res_calib_cosmic();
1260 test_cr2res_calib_flat();
1261 test_cr2res_calib_dark();
1262 test_cr2res_calib_bpm();
1263 test_cr2res_calib_detlin();
1264
1265 return cpl_test_end(0);
1266}
cr2res_calib_image
hdrl_image * cr2res_calib_image(const hdrl_image *in, int chip, int clean_bad, int subtract_nolight_rows, int subtract_interorder_column, int cosmics_corr, const cpl_frame *flat, const cpl_frame *dark, const cpl_frame *bpm, const cpl_frame *detlin, double dit, int ndit, int ndsamples)
The images calibration routine for a given chip.
Definition cr2res_calib.c:158
cr2res_io_save_MASTER_FLAT
int cr2res_io_save_MASTER_FLAT(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, hdrl_image **master_flats, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a MASTER_FLAT.
Definition cr2res_io.c:1669
cr2res_io_save_MASTER_DARK
int cr2res_io_save_MASTER_DARK(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, hdrl_image **master_darks, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a MASTER_DARK.
Definition cr2res_io.c:1521
cr2res_io_save_DETLIN_COEFFS
int cr2res_io_save_DETLIN_COEFFS(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, hdrl_imagelist **coeffs, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a DETLIN COEFFS.
Definition cr2res_io.c:1553
cr2res_io_save_BPM
int cr2res_io_save_BPM(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, cpl_image **bpms, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a BPM.
Definition cr2res_io.c:1585
hdrl_image_is_rejected
int hdrl_image_is_rejected(hdrl_image *self, cpl_size xpos, cpl_size ypos)
return if pixel is marked bad
Definition hdrl_image.c:445
hdrl_image_set_pixel
cpl_error_code hdrl_image_set_pixel(hdrl_image *self, cpl_size xpos, cpl_size ypos, hdrl_value value)
set pixel values of hdrl_image
Definition hdrl_image.c:594
hdrl_image_duplicate
hdrl_image * hdrl_image_duplicate(const hdrl_image *himg)
copy hdrl_image
Definition hdrl_image.c:391
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_add_scalar
cpl_error_code hdrl_image_add_scalar(hdrl_image *self, hdrl_value value)
Elementwise addition of a scalar to an image.
Definition hdrl_image_math.c:117
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_count_rejected
cpl_size hdrl_image_count_rejected(const hdrl_image *self)
return number of rejected pixels
Definition hdrl_image.c:461
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_reject
cpl_error_code hdrl_image_reject(hdrl_image *self, cpl_size xpos, cpl_size ypos)
mark pixel as bad
Definition hdrl_image.c:427
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_new
hdrl_image * hdrl_image_new(cpl_size nx, cpl_size ny)
create new zero filled hdrl image
Definition hdrl_image.c:311
hdrl_image_delete
void hdrl_image_delete(hdrl_image *himg)
delete hdrl_image
Definition hdrl_image.c:379
hdrl_imagelist_unset
hdrl_image * hdrl_imagelist_unset(hdrl_imagelist *himlist, cpl_size pos)
Remove an image from an imagelist.
Definition hdrl_imagelist_io.c:347
hdrl_imagelist_set
cpl_error_code hdrl_imagelist_set(hdrl_imagelist *himlist, hdrl_image *himg, cpl_size pos)
Insert an image into an imagelist.
Definition hdrl_imagelist_io.c:274
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_new
hdrl_imagelist * hdrl_imagelist_new(void)
Create an empty imagelist.
Definition hdrl_imagelist_io.c:84
Generated on Tue Oct 13 2026 12:27:25 for CR2RE Pipeline Reference Manual by doxygen 1.13.2