IIINSTRUMENT Pipeline Reference Manual 1.5.16
sofi_img_dark.c
1/* $Id: sofi_img_dark.c,v 1.27 2013-03-12 08:04:53 llundin Exp $
2 *
3 * This file is part of the SOFI Pipeline
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-12 08:04:53 $
24 * $Revision: 1.27 $
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 <cpl.h>
37#include <math.h>
38
39#include "irplib_utils.h"
40
41#include "sofi_utils.h"
42#include "sofi_pfits.h"
43#include "sofi_dfs.h"
44
45/*-----------------------------------------------------------------------------
46 Functions prototypes
47 -----------------------------------------------------------------------------*/
48
49static int sofi_img_dark_create(cpl_plugin *);
50static int sofi_img_dark_exec(cpl_plugin *);
51static int sofi_img_dark_destroy(cpl_plugin *);
52static int sofi_img_dark(cpl_parameterlist *, cpl_frameset *);
53
54static int sofi_img_dark_avg_reduce(cpl_frameset *, cpl_image **);
55static cpl_matrix * sofi_img_dark_ron_reduce(cpl_frameset *);
56static int sofi_img_dark_compare(const cpl_frame *, const cpl_frame *);
57static int sofi_img_dark_save(cpl_image *, cpl_matrix *, int, cpl_frameset *,
58 cpl_parameterlist *, cpl_frameset *);
59
60/*-----------------------------------------------------------------------------
61 Static variables
62 -----------------------------------------------------------------------------*/
63
64static struct {
65 /* Inputs */
66 int hsize;
67 int nsamples;
68 /* Outputs */
69 double dark_med;
70 double dark_stdev;
71} sofi_img_dark_config;
72
73static char sofi_img_dark_description[] =
74"sofi_img_dark -- SOFI imaging dark recipe.\n"
75"The files listed in the Set Of Frames (sof-file) must be tagged:\n"
76"raw-file.fits "SOFI_IMG_DARK_RAW"\n";
77
78/*-----------------------------------------------------------------------------
79 Functions code
80 -----------------------------------------------------------------------------*/
81
82/*----------------------------------------------------------------------------*/
90/*----------------------------------------------------------------------------*/
91int cpl_plugin_get_info(cpl_pluginlist * list)
92{
93 cpl_recipe * recipe = cpl_calloc(1, sizeof(*recipe));
94 cpl_plugin * plugin = &recipe->interface;
95
96 cpl_plugin_init(plugin,
97 CPL_PLUGIN_API,
98 SOFI_BINARY_VERSION,
99 CPL_PLUGIN_TYPE_RECIPE,
100 "sofi_img_dark",
101 "Dark recipe",
102 sofi_img_dark_description,
103 "Yves Jung",
104 "yjung@eso.org",
105 sofi_get_license(),
106 sofi_img_dark_create,
107 sofi_img_dark_exec,
108 sofi_img_dark_destroy);
109
110 cpl_pluginlist_append(list, plugin);
111
112 return 0;
113}
114
115/*----------------------------------------------------------------------------*/
124/*----------------------------------------------------------------------------*/
125static int sofi_img_dark_create(cpl_plugin * plugin)
126{
127 cpl_recipe * recipe;
128 cpl_parameter * p;
129
130 /* Get the recipe out of the plugin */
131 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
132 recipe = (cpl_recipe *)plugin;
133 else return CPL_ERROR_UNSPECIFIED;
134
135 /* Create the parameters list in the cpl_recipe object */
136 recipe->parameters = cpl_parameterlist_new();
137
138 /* Fill the parameters list */
139 /* --nsamples */
140 p = cpl_parameter_new_value("sofi.sofi_img_dark.nsamples",
141 CPL_TYPE_INT, "number of samples for RON computation",
142 "sofi.sofi_img_dark", 100);
143 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "nsamples");
144 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
145 cpl_parameterlist_append(recipe->parameters, p);
146 /* --hsize */
147 p = cpl_parameter_new_value("sofi.sofi_img_dark.hsize",
148 CPL_TYPE_INT, "half size of the window for RON computation",
149 "sofi.sofi_img_dark", 6);
150 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "hsize");
151 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
152 cpl_parameterlist_append(recipe->parameters, p);
153
154 /* Return */
155 return 0;
156}
157
158/*----------------------------------------------------------------------------*/
164/*----------------------------------------------------------------------------*/
165static int sofi_img_dark_exec(cpl_plugin * plugin)
166{
167 cpl_recipe * recipe;
168
169 /* Get the recipe out of the plugin */
170 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
171 recipe = (cpl_recipe *)plugin;
172 else return CPL_ERROR_UNSPECIFIED;
173
174 return sofi_img_dark(recipe->parameters, recipe->frames);
175}
176
177/*----------------------------------------------------------------------------*/
183/*----------------------------------------------------------------------------*/
184static int sofi_img_dark_destroy(cpl_plugin * plugin)
185{
186 cpl_recipe * recipe;
187
188 /* Get the recipe out of the plugin */
189 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
190 recipe = (cpl_recipe *)plugin;
191 else return CPL_ERROR_UNSPECIFIED;
192
193 cpl_parameterlist_delete(recipe->parameters);
194 return 0;
195}
196
197/*----------------------------------------------------------------------------*/
204/*----------------------------------------------------------------------------*/
205static int sofi_img_dark(
206 cpl_parameterlist * parlist,
207 cpl_frameset * framelist)
208{
209 cpl_parameter * par;
210 cpl_size nsets;
211 cpl_size * selection;
212 cpl_frameset * f_one;
213 cpl_image * avg;
214 cpl_matrix * rons;
215 int i;
216
217 /* Retrieve input parameters */
218 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_dark.hsize");
219 sofi_img_dark_config.hsize = cpl_parameter_get_int(par);
220 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_dark.nsamples");
221 sofi_img_dark_config.nsamples = cpl_parameter_get_int(par);
222
223 /* Identify the RAW and CALIB frames in the input frameset */
224 if (sofi_dfs_set_groups(framelist)) {
225 cpl_msg_error(cpl_func, "Cannot identify RAW and CALIB frames");
226 return CPL_ERROR_UNSPECIFIED;
227 }
228
229 /* Labelize all input frames */
230 cpl_msg_info(cpl_func, "Identify the different settings");
231 selection=cpl_frameset_labelise(framelist, sofi_img_dark_compare, &nsets);
232 if (selection == NULL) {
233 cpl_msg_error(cpl_func, "Cannot labelise input frames");
234 return CPL_ERROR_UNSPECIFIED;
235 }
236
237 /* Extract settings and reduce each of them */
238 for (i=0; i<nsets; i++) {
239 /* Initialise */
240 sofi_img_dark_config.dark_med = -1.0;
241 sofi_img_dark_config.dark_stdev = -1.0;
242 avg = NULL;
243
244 /* Reduce data set nb i */
245 cpl_msg_info(cpl_func, "Reduce data set no %d out of %d", i+1,
246 (int)nsets);
247 cpl_msg_indent_more();
248 f_one = cpl_frameset_extract(framelist, selection, i);
249
250 /* At least 2 frames required */
251 if (cpl_frameset_get_size(f_one) < 2) {
252 cpl_msg_warning(cpl_func, "Setting %d skipped (not enough frames)",
253 i+1);
254 cpl_frameset_delete(f_one);
255 cpl_free(selection);
256 return CPL_ERROR_UNSPECIFIED;
257 } else {
258 /* AVG part */
259 cpl_msg_info(cpl_func, "Compute the master dark");
260 cpl_msg_indent_more();
261 if (sofi_img_dark_avg_reduce(f_one, &avg)) {
262 cpl_msg_warning(cpl_func, "Cannot reduce set number %d", i+1);
263 }
264 cpl_msg_indent_less();
265 /* RON part */
266 cpl_msg_info(cpl_func, "Compute the read-out noise");
267 cpl_msg_indent_more();
268 if ((rons = sofi_img_dark_ron_reduce(f_one)) == NULL) {
269 cpl_msg_warning(cpl_func, "Cannot reduce set number %d", i+1);
270 }
271 cpl_msg_indent_less();
272 /* Save the products */
273 if (sofi_img_dark_save(avg, rons, i+1, f_one, parlist,
274 framelist) !=0 ){
275 cpl_msg_error(cpl_func, "Cannot save the products");
276 if (avg) cpl_image_delete(avg);
277 if (rons) cpl_matrix_delete(rons);
278 cpl_frameset_delete(f_one);
279 cpl_free(selection);
280 return CPL_ERROR_UNSPECIFIED;
281 }
282 if (rons) cpl_matrix_delete(rons);
283 }
284 if (avg) cpl_image_delete(avg);
285 cpl_frameset_delete(f_one);
286 cpl_msg_indent_less();
287 }
288
289 /* Free and return */
290 cpl_free(selection);
291 return 0;
292}
293
294/*----------------------------------------------------------------------------*/
301/*----------------------------------------------------------------------------*/
302static int sofi_img_dark_avg_reduce(
303 cpl_frameset * framelist,
304 cpl_image ** avg)
305{
306 cpl_imagelist * iset;
307 cpl_vector * medians;
308 cpl_image * image;
309 int i;
310
311 /* Test entries */
312 if (framelist == NULL) return CPL_ERROR_UNSPECIFIED;
313
314 /* Load the image set */
315 if ((iset = cpl_imagelist_load_frameset(framelist, CPL_TYPE_FLOAT, 1,
316 0)) == NULL) {
317 cpl_msg_error(cpl_func, "Cannot load the data");
318 return CPL_ERROR_UNSPECIFIED;
319 }
320
321 /* Average it to the master dark */
322 if ((*avg = cpl_imagelist_collapse_create(iset)) == NULL) {
323 cpl_msg_error(cpl_func, "Cannot average the data set");
324 cpl_imagelist_delete(iset);
325 return CPL_ERROR_UNSPECIFIED;
326 }
327
328 /* Compute mean-rms of the median values */
329 medians = cpl_vector_new(cpl_imagelist_get_size(iset));
330 for (i=0; i<cpl_imagelist_get_size(iset); i++) {
331 image = cpl_imagelist_get(iset, i);
332 cpl_vector_set(medians, i, cpl_image_get_median(image));
333 }
334 cpl_imagelist_delete(iset);
335 sofi_img_dark_config.dark_med = cpl_vector_get_mean(medians);
336 sofi_img_dark_config.dark_stdev = cpl_vector_get_stdev(medians);
337
338 /* Free and Return */
339 cpl_vector_delete(medians);
340 return 0;
341}
342
343/*----------------------------------------------------------------------------*/
354/*----------------------------------------------------------------------------*/
355static cpl_matrix * sofi_img_dark_ron_reduce(cpl_frameset * framelist)
356{
357 cpl_frame * cur_frame;
358 cpl_propertylist * plist;
359 cpl_imagelist * iset;
360 cpl_matrix * rons;
361 cpl_image * tmp_im;
362 double rms;
363 double norm;
364 int ndit;
365 cpl_size zone_def[4];
366 int i;
367
368 /* Test entries */
369 if (framelist == NULL) return NULL;
370
371 /* Load the current set */
372 if ((iset = cpl_imagelist_load_frameset(framelist, CPL_TYPE_FLOAT, 1,
373 0)) == NULL) {
374 cpl_msg_error(cpl_func, "Cannot load the data");
375 return NULL;
376 }
377
378 /* Create the matrix */
379 rons = cpl_matrix_new(cpl_imagelist_get_size(iset)-1, 4);
380
381 /* Loop on all pairs */
382 for (i=0; i<cpl_imagelist_get_size(iset)-1; i++) {
383 cpl_msg_info(cpl_func, "Pair number %d", i+1);
384 /* Get the norm factor */
385 if (cpl_error_get_code()) {
386 cpl_matrix_delete(rons);
387 cpl_imagelist_delete(iset);
388 return NULL;
389 }
390 cur_frame = cpl_frameset_get_position(framelist, i);
391 plist = cpl_propertylist_load(cpl_frame_get_filename(cur_frame),
392 0);
393 ndit = sofi_pfits_get_ndit(plist);
394 cpl_propertylist_delete(plist);
395 if (cpl_error_get_code()) {
396 cpl_msg_error(cpl_func, "Cannot get the NDIT");
397 cpl_matrix_delete(rons);
398 cpl_imagelist_delete(iset);
399 return NULL;
400 }
401 norm = 0.5 * ndit;
402 norm = sqrt(norm);
403
404 /* Compute the current subtracted image */
405 if ((tmp_im = cpl_image_subtract_create(
406 cpl_imagelist_get(iset, i),
407 cpl_imagelist_get(iset, i+1))) == NULL) {
408 cpl_msg_error(cpl_func, "Cannot subtract the images");
409 cpl_imagelist_delete(iset);
410 cpl_matrix_delete(rons);
411 return NULL;
412 }
413
414 /* Get measurement for lower-left quadrant */
415 zone_def[0] = 1;
416 zone_def[1] = cpl_image_get_size_x(tmp_im)/2;
417 zone_def[2] = 1;
418 zone_def[3] = cpl_image_get_size_y(tmp_im)/2;
419 cpl_flux_get_noise_window(tmp_im, zone_def,
420 sofi_img_dark_config.hsize,
421 sofi_img_dark_config.nsamples, &rms, NULL);
422 cpl_matrix_set(rons, i, 0, rms * norm);
423 cpl_msg_info(cpl_func, "RON in LL quadrant: %g", rms * norm);
424
425 /* Get measurement for lower-right quadrant */
426 zone_def[0] = cpl_image_get_size_x(tmp_im)/2 + 1;
427 zone_def[1] = cpl_image_get_size_x(tmp_im);
428 zone_def[2] = 1;
429 zone_def[3] = cpl_image_get_size_y(tmp_im)/2;
430 cpl_flux_get_noise_window(tmp_im, zone_def,
431 sofi_img_dark_config.hsize,
432 sofi_img_dark_config.nsamples, &rms, NULL);
433 cpl_matrix_set(rons, i, 1, rms * norm);
434 cpl_msg_info(cpl_func, "RON in LR quadrant: %g", rms * norm);
435
436 /* Get measurement for upper-left quadrant */
437 zone_def[0] = 1;
438 zone_def[1] = cpl_image_get_size_x(tmp_im)/2;
439 zone_def[2] = cpl_image_get_size_y(tmp_im)/2 + 1;
440 zone_def[3] = cpl_image_get_size_y(tmp_im);
441 cpl_flux_get_noise_window(tmp_im, zone_def,
442 sofi_img_dark_config.hsize,
443 sofi_img_dark_config.nsamples, &rms, NULL);
444 cpl_matrix_set(rons, i, 2, rms * norm);
445 cpl_msg_info(cpl_func, "RON in UL quadrant: %g", rms * norm);
446
447 /* Get measurement for upper-right quadrant */
448 zone_def[0] = cpl_image_get_size_x(tmp_im)/2 + 1;
449 zone_def[1] = cpl_image_get_size_x(tmp_im);
450 zone_def[2] = cpl_image_get_size_y(tmp_im)/2 + 1;
451 zone_def[3] = cpl_image_get_size_y(tmp_im);
452 cpl_flux_get_noise_window(tmp_im, zone_def,
453 sofi_img_dark_config.hsize,
454 sofi_img_dark_config.nsamples, &rms, NULL);
455 cpl_matrix_set(rons, i, 3, rms * norm);
456 cpl_msg_info(cpl_func, "RON in UR quadrant: %g", rms * norm);
457 cpl_image_delete(tmp_im);
458 }
459
460 /* Free and return */
461 cpl_imagelist_delete(iset);
462 return rons;
463}
464
465/*----------------------------------------------------------------------------*/
476/*----------------------------------------------------------------------------*/
477static int sofi_img_dark_save(
478 cpl_image * avg,
479 cpl_matrix * rons,
480 int set_nb,
481 cpl_frameset * set,
482 cpl_parameterlist * parlist,
483 cpl_frameset * set_tot)
484{
485 cpl_propertylist * plist;
486 cpl_propertylist * qclist;
487 cpl_propertylist * paflist;
488 const cpl_frame * ref_frame;
489 char qc_str[128];
490 char * filename;
491 int i;
492
493 /* Get the QC params in qclist */
494 qclist = cpl_propertylist_new();
495 cpl_propertylist_append_double(qclist, "ESO QC DARKMED",
496 sofi_img_dark_config.dark_med);
497 cpl_propertylist_append_double(qclist, "ESO QC DARKSTDEV",
498 sofi_img_dark_config.dark_stdev);
499 if (rons != NULL) {
500 for (i=0; i<cpl_matrix_get_nrow(rons); i++) {
501 sprintf(qc_str, "ESO QC LL RON%d", i+1);
502 cpl_propertylist_append_double(qclist, qc_str,
503 cpl_matrix_get(rons, i, 0));
504 sprintf(qc_str, "ESO QC LR RON%d", i+1);
505 cpl_propertylist_append_double(qclist, qc_str,
506 cpl_matrix_get(rons, i, 1));
507 sprintf(qc_str, "ESO QC UL RON%d", i+1);
508 cpl_propertylist_append_double(qclist, qc_str,
509 cpl_matrix_get(rons, i, 2));
510 sprintf(qc_str, "ESO QC UR RON%d", i+1);
511 cpl_propertylist_append_double(qclist, qc_str,
512 cpl_matrix_get(rons, i, 3));
513 }
514 }
515
516 /* Write the average image */
517 filename = cpl_sprintf("sofi_img_dark_set%02d_avg.fits", set_nb);
518 irplib_dfs_save_image(set_tot,
519 parlist,
520 set,
521 avg,
522 CPL_BPP_IEEE_FLOAT,
523 "sofi_img_dark",
524 SOFI_IMG_DARK_AVG,
525 qclist,
526 NULL,
527 PACKAGE "/" PACKAGE_VERSION,
528 filename);
529 cpl_free(filename);
530
531 /* Get the reference frame */
532 ref_frame = irplib_frameset_get_first_from_group(set, CPL_FRAME_GROUP_RAW);
533
534 /* Get FITS header from reference file */
535 if ((plist=cpl_propertylist_load(cpl_frame_get_filename(ref_frame),
536 0)) == NULL) {
537 cpl_msg_error(cpl_func, "getting header from reference frame");
538 cpl_propertylist_delete(qclist);
539 return CPL_ERROR_UNSPECIFIED;
540 }
541
542 /* Get the keywords for the paf file */
543 paflist = cpl_propertylist_new();
544 cpl_propertylist_copy_property_regexp(paflist, plist,
545 "^(ARCFILE|MJD-OBS|ESO TPL ID|ESO DPR TECH|DATE-OBS|ESO DET DIT|"
546 "ESO DET NDIT|ESO DET NCORRS|ESO DET NDSAMPLES|ESO DET MODE NAME)$", 0);
547 cpl_propertylist_delete(plist);
548
549 /* Copy the QC in paflist */
550 cpl_propertylist_copy_property_regexp(paflist, qclist, ".", 0);
551 cpl_propertylist_delete(qclist);
552
553 /* Save the PAF file */
554 filename = cpl_sprintf("sofi_img_dark_set%02d.paf", set_nb);
555 cpl_dfs_save_paf("SOFI",
556 "sofi_img_dark",
557 paflist,
558 filename);
559 cpl_free(filename);
560 cpl_propertylist_delete(paflist);
561 return 0;
562}
563
564/*----------------------------------------------------------------------------*/
571/*----------------------------------------------------------------------------*/
572static int sofi_img_dark_compare(
573 const cpl_frame * frame1,
574 const cpl_frame * frame2)
575{
576 int comparison;
577 cpl_propertylist * plist1;
578 cpl_propertylist * plist2;
579 double dval1, dval2;
580 int ival1, ival2;
581
582 /* Test entries */
583 if (frame1==NULL || frame2==NULL) return CPL_ERROR_UNSPECIFIED;
584
585 /* Get property lists */
586 if ((plist1=cpl_propertylist_load(cpl_frame_get_filename(frame1),
587 0)) == NULL) {
588 cpl_msg_error(cpl_func, "getting header from reference frame");
589 return CPL_ERROR_UNSPECIFIED;
590 }
591 if ((plist2=cpl_propertylist_load(cpl_frame_get_filename(frame2),
592 0)) == NULL) {
593 cpl_msg_error(cpl_func, "getting header from reference frame");
594 cpl_propertylist_delete(plist1);
595 return CPL_ERROR_UNSPECIFIED;
596 }
597
598 /* Test status */
599 if (cpl_error_get_code()) {
600 cpl_propertylist_delete(plist1);
601 cpl_propertylist_delete(plist2);
602 return CPL_ERROR_UNSPECIFIED;
603 }
604
605 /* Compare exposure time */
606 comparison = 1;
607 dval1 = sofi_pfits_get_dit(plist1);
608 dval2 = sofi_pfits_get_dit(plist2);
609 if (cpl_error_get_code()) {
610 cpl_msg_error(cpl_func, "cannot get exposure time");
611 cpl_propertylist_delete(plist1);
612 cpl_propertylist_delete(plist2);
613 return CPL_ERROR_UNSPECIFIED;
614 }
615 if (fabs(dval1-dval2) > 1e-5) comparison = 0;
616
617 /* Compare NDIT */
618 ival1 = sofi_pfits_get_ndit(plist1);
619 ival2 = sofi_pfits_get_ndit(plist2);
620 if (cpl_error_get_code()) {
621 cpl_msg_error(cpl_func, "cannot get NDIT");
622 cpl_propertylist_delete(plist1);
623 cpl_propertylist_delete(plist2);
624 return CPL_ERROR_UNSPECIFIED;
625 }
626 if (ival1 != ival2) comparison = 0;
627
628 /* Compare the readout mode */
629 ival1 = sofi_pfits_get_rom(plist1);
630 ival2 = sofi_pfits_get_rom(plist2);
631 if (cpl_error_get_code()) {
632 cpl_msg_error(cpl_func, "cannot get read-out mode");
633 cpl_propertylist_delete(plist1);
634 cpl_propertylist_delete(plist2);
635 return CPL_ERROR_UNSPECIFIED;
636 }
637 if (ival1 != ival2) comparison = 0;
638
639 /* Free and return */
640 cpl_propertylist_delete(plist1);
641 cpl_propertylist_delete(plist2);
642 return comparison;
643}
644
sofi_dfs_set_groups
int sofi_dfs_set_groups(cpl_frameset *set)
Set the group as RAW or CALIB in a frameset.
Definition sofi_dfs.c:60
sofi_pfits_get_dit
double sofi_pfits_get_dit(const cpl_propertylist *plist)
find out the DIT value
Definition sofi_pfits.c:195
sofi_pfits_get_rom
int sofi_pfits_get_rom(const cpl_propertylist *plist)
find out the read out mode
Definition sofi_pfits.c:486
sofi_pfits_get_ndit
int sofi_pfits_get_ndit(const cpl_propertylist *plist)
find out the NDIT keyword
Definition sofi_pfits.c:375
sofi_get_license
const char * sofi_get_license(void)
Get the pipeline copyright and license.
Definition sofi_utils.c:60
Generated on Thu Jun 18 2026 12:20:44 for IIINSTRUMENT Pipeline Reference Manual by doxygen 1.13.2