CR2RE Pipeline Reference Manual 1.7.0
irplib_calib.c
1/* $Id: irplib_calib.c,v 1.19 2013-03-01 10:26:22 llundin Exp $
2 *
3 * This file is part of the irplib package
4 * Copyright (C) 2002,2003 European Southern Observatory
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
19 */
20
21/*
22 * $Author: llundin $
23 * $Date: 2013-03-01 10:26:22 $
24 * $Revision: 1.19 $
25 * $Name: not supported by cvs2svn $
26 */
27
28#ifdef HAVE_CONFIG_H
29#include <config.h>
30#endif
31
32/*-----------------------------------------------------------------------------
33 Includes
34 -----------------------------------------------------------------------------*/
35
36#include "irplib_calib.h"
37
38#include <math.h>
39
40/*-----------------------------------------------------------------------------
41 Static Function Prototypes
42 -----------------------------------------------------------------------------*/
43
44static int irplib_get_clean_mean_window(cpl_image *img,
45 const int llx,
46 const int lly,
47 const int urx,
48 int ury,
49 const int kappa,
50 const int nclip,
51 double *clean_mean,
52 double *clean_stdev);
53
54static double irplib_pfits_get_dit(const cpl_propertylist *plist);
55static double irplib_pfits_get_exp_time(const cpl_propertylist *plist);
56/*----------------------------------------------------------------------------*/
60/*----------------------------------------------------------------------------*/
61
63/*---------------------------------------------------------------------------*/
72/*---------------------------------------------------------------------------*/
73static double
74irplib_pfits_get_dit(const cpl_propertylist *plist)
75{
76 cpl_errorstate prestate = cpl_errorstate_get();
77 double dit = cpl_propertylist_get_double(plist, "ESO DET DIT");
78
79 if (!cpl_errorstate_is_equal(prestate)) {
80 /* Key not present (or with wrong type, or something) */
81 cpl_errorstate prestate2 = cpl_errorstate_get();
82
83 dit = cpl_propertylist_get_double(plist, "ESO DET SEQ1 DIT");
84
85 if (cpl_errorstate_is_equal(prestate2)) {
86 /* Key present (with expected type): recover */
87 cpl_errorstate_set(prestate);
88 }
89 else {
90 cpl_error_set_where(cpl_func); /* Propagate error */
91 }
92 }
93
94 return dit;
95}
96
97/*---------------------------------------------------------------------------*/
105/*---------------------------------------------------------------------------*/
106static double
107irplib_pfits_get_exp_time(const cpl_propertylist *plist)
108{
109 return cpl_propertylist_get_double(plist, "EXPTIME");
110}
111
112
126
127
128static int
129irplib_get_clean_mean_window(cpl_image *img,
130 const int llx,
131 const int lly,
132 const int urx,
133 int ury,
134 const int kappa,
135 const int nclip,
136 double *clean_mean,
137 double *clean_stdev)
138{
139 double mean = 0;
140 double stdev = 0;
141 cpl_image *tmp = NULL;
142 cpl_stats *stats = NULL;
143 int i = 0;
144
145 tmp = cpl_image_extract(img, llx, lly, urx, ury);
146 cpl_image_accept_all(tmp);
147 for (i = 0; i < nclip; i++) {
148 double threshold = 0;
149 double lo_cut = 0;
150 double hi_cut = 0;
151 cpl_mask *mask = NULL;
152
153 cpl_stats_delete(stats);
154 stats = cpl_stats_new_from_image(tmp, CPL_STATS_MEAN | CPL_STATS_STDEV);
155 mean = cpl_stats_get_mean(stats);
156 stdev = cpl_stats_get_stdev(stats);
157
158 threshold = kappa * stdev;
159 lo_cut = mean - threshold;
160 hi_cut = mean + threshold;
161
162 cpl_image_accept_all(tmp);
163 mask = cpl_mask_threshold_image_create(tmp, lo_cut, hi_cut);
164
165 cpl_mask_not(mask);
166 cpl_image_reject_from_mask(tmp, mask);
167 cpl_mask_delete(mask);
168 }
169 *clean_mean = mean;
170 *clean_stdev = stdev;
171 cpl_image_delete(tmp);
172 cpl_stats_delete(stats);
173
174 return 0;
175}
176
177
178/*---------------------------------------------------------------------------*/
193/*---------------------------------------------------------------------------*/
194
195
196cpl_table *
197irplib_compute_gain(cpl_frameset *son,
198 cpl_frameset *sof,
199 int *zone,
200 const int kappa,
201 const int nclip)
202{
203 cpl_frame *frm = NULL;
204
205 cpl_table *res_tbl = NULL;
206 cpl_vector *dit_on = NULL;
207 cpl_vector *dit_of = NULL;
208 cpl_vector *exptime_on = NULL;
209 cpl_vector *exptime_of = NULL;
210
211 int non = 0;
212 int nof = 0;
213 int nfr = 0;
214 int llx;
215 int lly;
216 int urx;
217 int ury;
218
219 const char *name = NULL;
220 int i = 0;
221
222 double dit_ref = 0;
223 double exptime_ref = 0;
224
225
226 non = cpl_frameset_get_size(son);
227 nof = cpl_frameset_get_size(sof);
228 nfr = (non <= nof) ? non : nof;
229
230 dit_on = cpl_vector_new(nfr);
231 dit_of = cpl_vector_new(nfr);
232 exptime_on = cpl_vector_new(nfr);
233 exptime_of = cpl_vector_new(nfr);
234
235 for (i = 0; i < nfr; i++) {
236 cpl_propertylist *plist = NULL;
237
238 frm = cpl_frameset_get_position(son, i);
239 name = cpl_frame_get_filename(frm);
240 plist = cpl_propertylist_load(name, 0);
241 dit_ref = irplib_pfits_get_dit(plist);
242 exptime_ref = (double)irplib_pfits_get_exp_time(plist);
243 cpl_propertylist_delete(plist);
244 cpl_vector_set(dit_on, i, dit_ref);
245 cpl_vector_set(exptime_on, i, exptime_ref);
246
247 frm = cpl_frameset_get_position(sof, i);
248 name = cpl_frame_get_filename(frm);
249 plist = cpl_propertylist_load(name, 0);
250 dit_ref = irplib_pfits_get_dit(plist);
251 exptime_ref = (double)irplib_pfits_get_exp_time(plist);
252 cpl_propertylist_delete(plist);
253 cpl_vector_set(dit_of, i, dit_ref);
254 cpl_vector_set(exptime_of, i, exptime_ref);
255 }
256
257
258 llx = zone[0];
259 lly = zone[1];
260 urx = zone[2];
261 ury = zone[3];
262
263
264 res_tbl = cpl_table_new(nfr);
265 cpl_table_new_column(res_tbl, "adu", CPL_TYPE_DOUBLE);
266 cpl_table_new_column(res_tbl, "gain", CPL_TYPE_DOUBLE);
267
268 for (i = 0; i < nfr; i++) {
269 cpl_image *img_on1 = NULL;
270 cpl_image *img_of1 = NULL;
271 int m = 0;
272
273 frm = cpl_frameset_get_position(son, i);
274 name = cpl_frame_get_filename(frm);
275 img_on1 = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
276
277 frm = cpl_frameset_get_position(sof, i);
278 name = cpl_frame_get_filename(frm);
279 img_of1 = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
280
281
282 dit_ref = cpl_vector_get(dit_on, i);
283 exptime_ref = cpl_vector_get(exptime_on, i);
284
285
286 for (m = 0; m < nfr; m++) {
287 if (m != i) {
288 double dit_tmp = 0;
289 double exptime_tmp = 0;
290
291 frm = cpl_frameset_get_position(son, m);
292 name = cpl_frame_get_filename(frm);
293 dit_tmp = cpl_vector_get(dit_on, m);
294 exptime_tmp = cpl_vector_get(exptime_on, m);
295 if (dit_tmp == dit_ref && exptime_tmp == exptime_ref) {
296 cpl_image *img_on2 = NULL;
297 cpl_image *img_on_dif = NULL;
298
299 cpl_image *img_of2 = NULL;
300 cpl_image *img_of_dif = NULL;
301
302 double avg_on1 = 0;
303 double avg_on2 = 0;
304 double avg_of1 = 0;
305 double avg_of2 = 0;
306 double avg_on_dif = 0;
307 double avg_of_dif = 0;
308 double std = 0;
309 double sig_on_dif = 0;
310 double sig_of_dif = 0;
311 double gain = 0;
312
313 img_on2 = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
314 frm = cpl_frameset_get_position(sof, m);
315 name = cpl_frame_get_filename(frm);
316 img_of2 = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
317
318 img_on_dif = cpl_image_subtract_create(img_on1, img_on2);
319 img_of_dif = cpl_image_subtract_create(img_of1, img_of2);
320
321 irplib_get_clean_mean_window(img_on1, llx, lly, urx, ury,
322 kappa, nclip, &avg_on1, &std);
323 irplib_get_clean_mean_window(img_on2, llx, lly, urx, ury,
324 kappa, nclip, &avg_on2, &std);
325 irplib_get_clean_mean_window(img_of1, llx, lly, urx, ury,
326 kappa, nclip, &avg_of1, &std);
327 irplib_get_clean_mean_window(img_of2, llx, lly, urx, ury,
328 kappa, nclip, &avg_of2, &std);
329 irplib_get_clean_mean_window(img_on_dif, llx, lly, urx, ury,
330 kappa, nclip, &avg_on_dif,
331 &sig_on_dif);
332 irplib_get_clean_mean_window(img_of_dif, llx, lly, urx, ury,
333 kappa, nclip, &avg_of_dif,
334 &sig_of_dif);
335
336 cpl_image_delete(img_on2);
337 cpl_image_delete(img_of2);
338 cpl_image_delete(img_on_dif);
339 cpl_image_delete(img_of_dif);
340
341 gain =
342 ((avg_on1 + avg_on2) - (avg_of1 + avg_of2)) /
343 ((sig_on_dif * sig_on_dif) - (sig_of_dif * sig_of_dif));
344
345 cpl_table_set_double(res_tbl, "gain", m, gain);
346 cpl_table_set_double(res_tbl, "adu", m,
347 ((avg_on1 + avg_on2) / 2 -
348 (avg_of1 + avg_of2) / 2));
349 }
350 }
351 }
352 cpl_image_delete(img_on1);
353 cpl_image_delete(img_of1);
354 }
355
356
357 cpl_vector_delete(dit_on);
358 cpl_vector_delete(dit_of);
359 cpl_vector_delete(exptime_on);
360 cpl_vector_delete(exptime_of);
361
362 return res_tbl;
363}
364
365/* --------------------------------------------------------------------------*/
375/*---------------------------------------------------------------------------*/
376
377
378cpl_table *
379irplib_compute_linearity(cpl_frameset *son, cpl_frameset *sof)
380{
381 int non = 0;
382 int nof = 0;
383 int nfr = 0;
384 int i = 0;
385 double med_dit = 0;
386 /*double avg_dit=0;*/
387
388 cpl_vector *vec_adl = NULL;
389 cpl_vector *vec_dit = NULL;
390 cpl_vector *vec_avg = NULL;
391 cpl_vector *vec_med = NULL;
392 cpl_vector *vec_avg_dit = NULL;
393 cpl_vector *vec_med_dit = NULL;
394
395 double dit = 0;
396 cpl_table *lin_tbl = NULL;
397
398
399 non = cpl_frameset_get_size(son);
400 nof = cpl_frameset_get_size(sof);
401 nfr = (non <= nof) ? non : nof;
402
403 lin_tbl = cpl_table_new(nfr);
404 cpl_table_new_column(lin_tbl, "med", CPL_TYPE_DOUBLE);
405 cpl_table_new_column(lin_tbl, "avg", CPL_TYPE_DOUBLE);
406 cpl_table_new_column(lin_tbl, "med_dit", CPL_TYPE_DOUBLE);
407 cpl_table_new_column(lin_tbl, "avg_dit", CPL_TYPE_DOUBLE);
408 cpl_table_new_column(lin_tbl, "dit", CPL_TYPE_DOUBLE);
409 vec_med = cpl_vector_new(nfr);
410 vec_avg = cpl_vector_new(nfr);
411 vec_med_dit = cpl_vector_new(nfr);
412 vec_avg_dit = cpl_vector_new(nfr);
413 vec_dit = cpl_vector_new(nfr);
414 vec_adl = cpl_vector_new(nfr);
415 for (i = 0; i < nfr; i++) {
416 cpl_frame *frm = NULL;
417
418 double med_on = 0;
419 double avg_on = 0;
420 double med_of = 0;
421 double avg_of = 0;
422
423 double med = 0;
424 double avg = 0;
425 double avg_dit = 0;
426
427 const char *name = NULL;
428 cpl_image *img = NULL;
429 cpl_propertylist *plist = NULL;
430
431 frm = cpl_frameset_get_position(son, i);
432 name = cpl_frame_get_filename(frm);
433 img = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
434 med_on = cpl_image_get_median(img);
435 avg_on = cpl_image_get_mean(img);
436 cpl_image_delete(img);
437
438 frm = cpl_frameset_get_position(sof, i);
439 name = cpl_frame_get_filename(frm);
440 img = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
441 med_of = cpl_image_get_median(img);
442 avg_of = cpl_image_get_mean(img);
443 cpl_image_delete(img);
444 med = med_on - med_of;
445 avg = avg_on - avg_of;
446 plist = cpl_propertylist_load(name, 0);
447 dit = (double)irplib_pfits_get_dit(plist);
448 cpl_propertylist_delete(plist);
449 avg_dit = avg / dit;
450 med_dit = med / dit;
451
452 cpl_vector_set(vec_dit, i, dit);
453 cpl_vector_set(vec_avg, i, avg);
454 cpl_vector_set(vec_med, i, med);
455 cpl_vector_set(vec_avg_dit, i, avg_dit);
456 cpl_vector_set(vec_med_dit, i, med_dit);
457
458
459 cpl_table_set_double(lin_tbl, "dit", i, dit);
460 cpl_table_set_double(lin_tbl, "med", i, med);
461 cpl_table_set_double(lin_tbl, "avg", i, avg);
462 cpl_table_set_double(lin_tbl, "med_dit", i, med_dit);
463 cpl_table_set_double(lin_tbl, "avg_dit", i, avg_dit);
464 }
465 cpl_table_new_column(lin_tbl, "adl", CPL_TYPE_DOUBLE);
466 med_dit = cpl_vector_get_mean(vec_med_dit);
467 /*avg_dit=cpl_vector_get_mean(vec_avg_dit);*/
468
469 for (i = 0; i < nfr; i++) {
470 int *status = 0;
471 dit = cpl_table_get_double(lin_tbl, "dit", i, status);
472 cpl_vector_set(vec_adl, i, dit * med_dit);
473 cpl_table_set_double(lin_tbl, "adl", i, dit * med_dit);
474 }
475
476
477 cpl_vector_delete(vec_dit);
478 cpl_vector_delete(vec_adl);
479 cpl_vector_delete(vec_avg);
480 cpl_vector_delete(vec_med);
481 cpl_vector_delete(vec_avg_dit);
482 cpl_vector_delete(vec_med_dit);
483
484
485 return lin_tbl;
486}
487
488
489/*----------------------------------------------------------------------------*/
498/*----------------------------------------------------------------------------*/
499int
500irplib_detlin_correct(cpl_imagelist *ilist,
501 const char *detlin_a,
502 const char *detlin_b,
503 const char *detlin_c)
504{
505 cpl_image *ima;
506 cpl_image *imb;
507 cpl_image *imc;
508 float *pima;
509 float *pimb;
510 float *pimc;
511 float *pdata;
512 int nx, ny, ni;
513 double coeff_1, coeff_2, val;
514 int i, j;
515
516 /* Test entries */
517 if (!ilist || !detlin_a || !detlin_b || !detlin_c)
518 return -1;
519
520 /* Load the 3 coeffs images */
521 ima = cpl_image_load(detlin_a, CPL_TYPE_FLOAT, 0, 0);
522 imb = cpl_image_load(detlin_b, CPL_TYPE_FLOAT, 0, 0);
523 imc = cpl_image_load(detlin_c, CPL_TYPE_FLOAT, 0, 0);
524 if (!ima || !imb || !imc) {
525 cpl_msg_error(cpl_func, "Cannot load the detlin images");
526 if (ima)
527 cpl_image_delete(ima);
528 if (imb)
529 cpl_image_delete(imb);
530 if (imc)
531 cpl_image_delete(imc);
532 return -1;
533 }
534 pima = cpl_image_get_data_float(ima);
535 pimb = cpl_image_get_data_float(imb);
536 pimc = cpl_image_get_data_float(imc);
537
538 /* Test sizes */
539 nx = cpl_image_get_size_x(cpl_imagelist_get(ilist, 0));
540 ny = cpl_image_get_size_y(cpl_imagelist_get(ilist, 0));
541 ni = cpl_imagelist_get_size(ilist);
542 if ((cpl_image_get_size_x(ima) != nx) ||
543 (cpl_image_get_size_x(imb) != nx) ||
544 (cpl_image_get_size_x(imc) != nx) ||
545 (cpl_image_get_size_y(ima) != ny) ||
546 (cpl_image_get_size_y(imb) != ny) ||
547 (cpl_image_get_size_y(imc) != ny)) {
548 cpl_msg_error(cpl_func, "Incompatible sizes");
549 cpl_image_delete(ima);
550 cpl_image_delete(imb);
551 cpl_image_delete(imc);
552 return -1;
553 }
554
555 /* Loop on pixels */
556 for (i = 0; i < nx * ny; i++) {
557 /* Compute the coefficients */
558 if (fabs(pima[i]) < 1e-30) {
559 coeff_1 = coeff_2 = (double)0.0;
560 }
561 else {
562 coeff_1 = (double)pimb[i] / (double)pima[i];
563 coeff_2 = (double)pimc[i] / (double)pima[i];
564 }
565 /* Correct this pixel in each plane */
566 for (j = 0; j < ni; j++) {
567 pdata = cpl_image_get_data_float(cpl_imagelist_get(ilist, j));
568 val = (double)pdata[i];
569 pdata[i] =
570 (float)(val + coeff_1 * val * val + coeff_2 * val * val * val);
571 }
572 }
573 /* Free and return */
574 cpl_image_delete(ima);
575 cpl_image_delete(imb);
576 cpl_image_delete(imc);
577 return 0;
578}
579
580/*----------------------------------------------------------------------------*/
589/*----------------------------------------------------------------------------*/
590int
591irplib_flat_dark_bpm_calib(cpl_imagelist *ilist,
592 const char *flat,
593 const char *dark,
594 const char *bpm)
595{
596 /* Test entries */
597 if (ilist == NULL)
598 return -1;
599
600 /* Dark correction */
601 if (dark != NULL) {
602 cpl_image *dark_image;
603 cpl_msg_info(cpl_func, "Subtract the dark to the images");
604 /* Load the dark image */
605 if ((dark_image = cpl_image_load(dark, CPL_TYPE_FLOAT, 0, 0)) == NULL) {
606 cpl_msg_error(cpl_func, "Cannot load the dark %s", dark);
607 return -1;
608 }
609 /* Apply the dark correction to the images */
610 if (cpl_imagelist_subtract_image(ilist, dark_image) != CPL_ERROR_NONE) {
611 cpl_msg_error(cpl_func, "Cannot apply the dark to the images");
612 cpl_image_delete(dark_image);
613 return -1;
614 }
615 cpl_image_delete(dark_image);
616 }
617
618 /* Flat-field correction */
619 if (flat != NULL) {
620 cpl_image *flat_image;
621 cpl_msg_info(cpl_func, "Divide the images by the flatfield");
622 /* Load the flat image */
623 if ((flat_image = cpl_image_load(flat, CPL_TYPE_FLOAT, 0, 0)) == NULL) {
624 cpl_msg_error(cpl_func, "Cannot load the flat field %s", flat);
625 return -1;
626 }
627 /* Apply the flatfield correction to the images */
628 if (cpl_imagelist_divide_image(ilist, flat_image) != CPL_ERROR_NONE) {
629 cpl_msg_error(cpl_func, "Cannot apply the flatfield to the images");
630 cpl_image_delete(flat_image);
631 return -1;
632 }
633 cpl_image_delete(flat_image);
634 }
635
636 /* Correct the bad pixels if requested */
637 if (bpm != NULL) {
638 cpl_mask *bpm_im_bin;
639 cpl_image *bpm_im_int;
640 int i;
641 cpl_msg_info(cpl_func, "Correct the bad pixels in the images");
642 /* Load the bad pixels image */
643 if ((bpm_im_int = cpl_image_load(bpm, CPL_TYPE_INT, 0, 0)) == NULL) {
644 cpl_msg_error(cpl_func, "Cannot load the bad pixel map %s", bpm);
645 return -1;
646 }
647 /* Convert the map from integer to binary */
648 bpm_im_bin = cpl_mask_threshold_image_create(bpm_im_int, -0.5, 0.5);
649 cpl_mask_not(bpm_im_bin);
650 cpl_image_delete(bpm_im_int);
651 /* Apply the bad pixels cleaning */
652 for (i = 0; i < cpl_imagelist_get_size(ilist); i++) {
653 cpl_image_reject_from_mask(cpl_imagelist_get(ilist, i), bpm_im_bin);
654 if (cpl_detector_interpolate_rejected(
655 cpl_imagelist_get(ilist, i)) != CPL_ERROR_NONE) {
656 cpl_msg_error(cpl_func, "Cannot clean the bad pixels in obj %d",
657 i + 1);
658 cpl_mask_delete(bpm_im_bin);
659 return -1;
660 }
661 }
662 cpl_mask_delete(bpm_im_bin);
663 }
664
665 /* Return */
666 return 0;
667}
668
irplib_compute_gain
cpl_table * irplib_compute_gain(cpl_frameset *son, cpl_frameset *sof, int *zone, const int kappa, const int nclip)
Computes the detector's gain.
Definition irplib_calib.c:197
irplib_compute_linearity
cpl_table * irplib_compute_linearity(cpl_frameset *son, cpl_frameset *sof)
Computes the detector's linearity.
Definition irplib_calib.c:379
irplib_flat_dark_bpm_calib
int irplib_flat_dark_bpm_calib(cpl_imagelist *ilist, const char *flat, const char *dark, const char *bpm)
Apply the calibration to the frames.
Definition irplib_calib.c:591
irplib_detlin_correct
int irplib_detlin_correct(cpl_imagelist *ilist, const char *detlin_a, const char *detlin_b, const char *detlin_c)
Apply the detector linearity correction.
Definition irplib_calib.c:500
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2