IIINSTRUMENT Pipeline Reference Manual 1.5.16
sofi_img_detlin.c
1/* $Id: sofi_img_detlin.c,v 1.16 2013-03-12 08:04:51 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:51 $
24 * $Revision: 1.16 $
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 <math.h>
37#include <cpl.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_detlin_create(cpl_plugin *) ;
50static int sofi_img_detlin_exec(cpl_plugin *) ;
51static int sofi_img_detlin_destroy(cpl_plugin *) ;
52static int sofi_img_detlin(cpl_parameterlist *, cpl_frameset *) ;
53
54static cpl_imagelist * sofi_img_detlin_load(cpl_frameset *, cpl_frameset *,
55 cpl_vector **) ;
56static int sofi_img_detlin_save(cpl_imagelist *, cpl_parameterlist *,
57 cpl_frameset *) ;
58
59/*-----------------------------------------------------------------------------
60 Static variables
61 -----------------------------------------------------------------------------*/
62
63static struct {
64 /* Inputs */
65 int force_flag ;
66
67 /* Outputs */
68 double lamp_stability ;
69} sofi_img_detlin_config ;
70
71static char sofi_img_detlin_description[] =
72"sofi_img_detlin -- SOFI imaging detector linearity recipe.\n"
73"The files listed in the Set Of Frames (sof-file) must be tagged:\n"
74"raw-file.fits "SOFI_IMG_DETLIN_LAMP_RAW" or\n"
75"raw-file.fits "SOFI_IMG_DETLIN_DARK_RAW"\n";
76
77/*-----------------------------------------------------------------------------
78 Functions code
79 -----------------------------------------------------------------------------*/
80
81/*----------------------------------------------------------------------------*/
89/*----------------------------------------------------------------------------*/
90int cpl_plugin_get_info(cpl_pluginlist * list)
91{
92 cpl_recipe * recipe = cpl_calloc(1, sizeof(*recipe)) ;
93 cpl_plugin * plugin = &recipe->interface ;
94
95 cpl_plugin_init(plugin,
96 CPL_PLUGIN_API,
97 SOFI_BINARY_VERSION,
98 CPL_PLUGIN_TYPE_RECIPE,
99 "sofi_img_detlin",
100 "Detector linearity recipe",
101 sofi_img_detlin_description,
102 "Yves Jung",
103 "yjung@eso.org",
105 sofi_img_detlin_create,
106 sofi_img_detlin_exec,
107 sofi_img_detlin_destroy) ;
108
109 cpl_pluginlist_append(list, plugin) ;
110
111 return 0;
112}
113
114/*----------------------------------------------------------------------------*/
123/*----------------------------------------------------------------------------*/
124static int sofi_img_detlin_create(cpl_plugin * plugin)
125{
126 cpl_recipe * recipe ;
127 cpl_parameter * p ;
128
129 /* Get the recipe out of the plugin */
130 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
131 recipe = (cpl_recipe *)plugin ;
132 else return CPL_ERROR_UNSPECIFIED;
133
134 /* Create the parameters list in the cpl_recipe object */
135 recipe->parameters = cpl_parameterlist_new() ;
136
137 /* Fill the parameters list */
138 /* --force */
139 p = cpl_parameter_new_value("sofi.sofi_img_detlin.force", CPL_TYPE_BOOL,
140 "flag to force th computation", "sofi.sofi_img_detlin",
141 FALSE) ;
142 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "force") ;
143 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV) ;
144 cpl_parameterlist_append(recipe->parameters, p) ;
145
146 /* Return */
147 return 0;
148}
149
150/*----------------------------------------------------------------------------*/
156/*----------------------------------------------------------------------------*/
157static int sofi_img_detlin_exec(cpl_plugin * plugin)
158{
159 cpl_recipe * recipe ;
160
161 /* Get the recipe out of the plugin */
162 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
163 recipe = (cpl_recipe *)plugin ;
164 else return CPL_ERROR_UNSPECIFIED;
165
166 return sofi_img_detlin(recipe->parameters, recipe->frames) ;
167}
168
169/*----------------------------------------------------------------------------*/
175/*----------------------------------------------------------------------------*/
176static int sofi_img_detlin_destroy(cpl_plugin * plugin)
177{
178 cpl_recipe * recipe ;
179
180 /* Get the recipe out of the plugin */
181 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
182 recipe = (cpl_recipe *)plugin ;
183 else return CPL_ERROR_UNSPECIFIED;
184
185 cpl_parameterlist_delete(recipe->parameters) ;
186 return 0 ;
187}
188
189/*----------------------------------------------------------------------------*/
196/*----------------------------------------------------------------------------*/
197static int sofi_img_detlin(
198 cpl_parameterlist * parlist,
199 cpl_frameset * framelist)
200{
201 cpl_parameter * par ;
202 cpl_frameset * darkframes ;
203 cpl_frameset * lampframes ;
204 cpl_imagelist * iset ;
205 cpl_vector * ditval ;
206 cpl_imagelist * fitres ;
207 cpl_image * error_im ;
208 int degree ;
209
210 /* Initialise */
211 par = NULL ;
212 degree = 2 ;
213
214 /* Retrieve input parameters */
215 /* Force flag */
216 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_detlin.force") ;
217 sofi_img_detlin_config.force_flag = cpl_parameter_get_bool(par) ;
218
219 /* Identify the RAW and CALIB frames in the input frameset */
220 if (sofi_dfs_set_groups(framelist)) {
221 cpl_msg_error(cpl_func, "Cannot identify RAW and CALIB frames") ;
222 return CPL_ERROR_UNSPECIFIED;
223 }
224
225 /* Retrieve raw frames */
226 if ((lampframes = sofi_extract_frameset(framelist,
227 SOFI_IMG_DETLIN_LAMP_RAW)) == NULL) {
228 cpl_msg_error(cpl_func, "Cannot find any lamp frame in input");
229 return CPL_ERROR_UNSPECIFIED;
230 }
231 if ((darkframes = sofi_extract_frameset(framelist,
232 SOFI_IMG_DETLIN_DARK_RAW)) == NULL) {
233 cpl_msg_error(cpl_func, "Cannot find any dark frame in input");
234 cpl_frameset_delete(lampframes) ;
235 return CPL_ERROR_UNSPECIFIED;
236 }
237
238 /* Load the data and check the DIT consistency */
239 cpl_msg_info(cpl_func, "Load the data") ;
240 cpl_msg_indent_more() ;
241 if ((iset = sofi_img_detlin_load(lampframes, darkframes, &ditval))==NULL) {
242 cpl_msg_error(cpl_func, "Cannot load the data") ;
243 cpl_frameset_delete(lampframes) ;
244 cpl_frameset_delete(darkframes) ;
245 cpl_msg_indent_less() ;
246 return CPL_ERROR_UNSPECIFIED;
247 }
248 cpl_frameset_delete(lampframes) ;
249 cpl_frameset_delete(darkframes) ;
250 cpl_msg_indent_less() ;
251
252 /* Create the error image */
253 error_im = cpl_image_duplicate(cpl_imagelist_get(iset, 0)) ;
254
255 /* Compute the linearity coefficients */
256 cpl_msg_info(cpl_func, "Compute the linearity coefficients") ;
257 cpl_msg_indent_more() ;
258 if ((fitres = cpl_fit_imagelist_polynomial(ditval, iset, 0, degree,
259 CPL_FALSE, CPL_TYPE_FLOAT, error_im)) == NULL) {
260 cpl_msg_error(cpl_func, "Cannot compute the linearity coefficients") ;
261 cpl_imagelist_delete(iset) ;
262 cpl_vector_delete(ditval) ;
263 cpl_image_delete(error_im) ;
264 cpl_msg_indent_less() ;
265 return CPL_ERROR_UNSPECIFIED;
266 }
267 cpl_msg_indent_less() ;
268 cpl_vector_delete(ditval) ;
269 cpl_imagelist_delete(iset) ;
270
271 /* Add the error image to the coefficients */
272 cpl_imagelist_set(fitres, error_im, degree+1) ;
273
274 /* Save the products */
275 cpl_msg_info(cpl_func, "Save the products") ;
276 cpl_msg_indent_more() ;
277 if (sofi_img_detlin_save(fitres, parlist, framelist)==-1) {
278 cpl_msg_error(cpl_func, "Cannot save the products") ;
279 cpl_imagelist_delete(fitres) ;
280 cpl_msg_indent_less() ;
281 return CPL_ERROR_UNSPECIFIED;
282 }
283 cpl_msg_indent_less() ;
284
285 /* Free and return */
286 cpl_imagelist_delete(fitres) ;
287 return 0 ;
288}
289
290/*----------------------------------------------------------------------------*/
298/*----------------------------------------------------------------------------*/
299static cpl_imagelist * sofi_img_detlin_load(
300 cpl_frameset * lamps,
301 cpl_frameset * darks,
302 cpl_vector ** ditvals)
303{
304 int nb_lamps ;
305 cpl_vector * selection ;
306 cpl_frame * frame ;
307 cpl_propertylist * propertylist ;
308 double dit_lamp, dit_dark ;
309 int dit_stab ;
310 cpl_imagelist * lamps_data ;
311 cpl_imagelist * darks_data ;
312 double * stab_levels ;
313 cpl_vector * dit_purged ;
314 double * pditvals ;
315 double * pdit_purged ;
316 int i, j ;
317
318 /* Check that there are as many lamp as darks */
319 if ((nb_lamps = cpl_frameset_get_size(lamps)) < 3) return NULL ;
320 if (cpl_frameset_get_size(darks) != nb_lamps) return NULL ;
321
322 /* Check out that they have consistent integration times */
323 cpl_msg_info(cpl_func, "Checking DIT consistency") ;
324 selection = cpl_vector_new(nb_lamps) ;
325 *ditvals = cpl_vector_new(nb_lamps) ;
326 pditvals = cpl_vector_get_data(*ditvals) ;
327 dit_stab = 0 ;
328 for (i=0 ; i<nb_lamps ; i++) {
329 /* Check if ok */
330 if (cpl_error_get_code()) {
331 cpl_vector_delete(selection) ;
332 cpl_vector_delete(*ditvals) ;
333 return NULL ;
334 }
335 /* DIT from LAMP */
336 frame = cpl_frameset_get_position(lamps, i) ;
337 propertylist=cpl_propertylist_load(cpl_frame_get_filename(frame), 0) ;
338 dit_lamp = (double)sofi_pfits_get_dit(propertylist) ;
339 cpl_propertylist_delete(propertylist) ;
340 if (cpl_error_get_code()) {
341 cpl_msg_error(cpl_func, "Cannot get DIT") ;
342 cpl_vector_delete(selection) ;
343 cpl_vector_delete(*ditvals) ;
344 return NULL ;
345 }
346 /* DIT from DARK */
347 frame = cpl_frameset_get_position(darks, i) ;
348 propertylist=cpl_propertylist_load(cpl_frame_get_filename(frame), 0) ;
349 dit_dark = (double)sofi_pfits_get_dit(propertylist) ;
350 cpl_propertylist_delete(propertylist) ;
351 if (cpl_error_get_code()) {
352 cpl_msg_error(cpl_func, "Cannot get DIT") ;
353 cpl_vector_delete(selection) ;
354 cpl_vector_delete(*ditvals) ;
355 return NULL ;
356 }
357 /* Check consistency */
358 if (fabs(dit_dark-dit_lamp) > 1e-3) {
359 cpl_msg_error(cpl_func, "DIT not consistent between LAMP and DARK");
360 cpl_vector_delete(selection) ;
361 cpl_vector_delete(*ditvals) ;
362 return NULL ;
363 }
364 pditvals[i] = dit_lamp ;
365 /* Set selection */
366 if (i==0) {
367 cpl_vector_set(selection, i, -1.0) ;
368 dit_stab ++ ;
369 } else {
370 if (fabs(dit_lamp - pditvals[0]) < 1e-5) {
371 cpl_vector_set(selection, i, -1.0) ;
372 dit_stab ++ ;
373 } else {
374 cpl_vector_set(selection, i, 1.0) ;
375 }
376 }
377 }
378
379 /* Check if there are enough DITs for stability check */
380 if (dit_stab < 2) {
381 cpl_msg_error(cpl_func, "Not enough frames for stability check") ;
382 cpl_vector_delete(selection) ;
383 cpl_vector_delete(*ditvals) ;
384 return NULL ;
385 }
386
387 /* Load the data and compute lamp-dark */
388 cpl_msg_info(cpl_func, "Compute the differences lamp - dark") ;
389 lamps_data = cpl_imagelist_load_frameset(lamps, CPL_TYPE_FLOAT, 1, 0) ;
390 darks_data = cpl_imagelist_load_frameset(darks, CPL_TYPE_FLOAT, 1, 0) ;
391 if (cpl_imagelist_subtract(lamps_data,darks_data) != CPL_ERROR_NONE) {
392 cpl_msg_error(cpl_func, "Cannot subtract the 2 image lists") ;
393 cpl_vector_delete(selection) ;
394 cpl_vector_delete(*ditvals) ;
395 if (lamps_data) cpl_imagelist_delete(lamps_data) ;
396 if (darks_data) cpl_imagelist_delete(darks_data) ;
397 return NULL ;
398 }
399 cpl_imagelist_delete(darks_data) ;
400
401 /* Check the lamp stability */
402 cpl_msg_info(cpl_func, "Check the lamp stability") ;
403 stab_levels = cpl_malloc(dit_stab * sizeof(double)) ;
404 j = 0 ;
405 for (i=0 ; i<nb_lamps ; i++) {
406 if (cpl_vector_get(selection, i) < 0) {
407 stab_levels[j] =
408 cpl_image_get_mean(cpl_imagelist_get(lamps_data, i)) ;
409 j++ ;
410 }
411 }
412
413 /* Compute the lamp stability */
414 sofi_img_detlin_config.lamp_stability = 0.0 ;
415 for (i=1 ; i<dit_stab ; i++) {
416 if ((fabs(stab_levels[i]-stab_levels[0]) / stab_levels[0]) >
417 sofi_img_detlin_config.lamp_stability)
418 sofi_img_detlin_config.lamp_stability =
419 fabs(stab_levels[i]-stab_levels[0]) / stab_levels[0] ;
420 }
421 cpl_free(stab_levels) ;
422
423 /* Check the lamp stability */
424 if (sofi_img_detlin_config.lamp_stability > 0.01) {
425 if (sofi_img_detlin_config.force_flag == 1) {
426 cpl_msg_warning(cpl_func,
427 "level difference #%d too high - proceed anyway",i+1);
428 } else {
429 cpl_msg_error(cpl_func, "level difference #%d too high", i+1);
430 cpl_vector_delete(selection) ;
431 cpl_vector_delete(*ditvals) ;
432 cpl_imagelist_delete(lamps_data) ;
433 return NULL ;
434 }
435 }
436
437 /* Purge the stability check data and DITs */
438 if (cpl_imagelist_erase(lamps_data, selection) != CPL_ERROR_NONE) {
439 cpl_msg_error(cpl_func, "cannot discard stability frames") ;
440 cpl_vector_delete(selection) ;
441 cpl_vector_delete(*ditvals) ;
442 cpl_imagelist_delete(lamps_data) ;
443 return NULL ;
444 }
445 dit_purged = cpl_vector_new(cpl_imagelist_get_size(lamps_data)) ;
446 pdit_purged = cpl_vector_get_data(dit_purged) ;
447 j = 0 ;
448 for (i=0 ; i<nb_lamps ; i++) {
449 if (cpl_vector_get(selection, i) > 0) {
450 pdit_purged[j] = pditvals[i] ;
451 j++ ;
452 }
453 }
454 cpl_vector_delete(*ditvals) ;
455 *ditvals = dit_purged ;
456
457 /* Free and eturn */
458 cpl_vector_delete(selection) ;
459 return lamps_data ;
460}
461
462/*----------------------------------------------------------------------------*/
470/*----------------------------------------------------------------------------*/
471static int sofi_img_detlin_save(
472 cpl_imagelist * fitres,
473 cpl_parameterlist * parlist,
474 cpl_frameset * set)
475{
476 const cpl_frame * ref_frame;
477 cpl_propertylist * plist ;
478 cpl_propertylist * qclist ;
479 cpl_propertylist * paflist ;
480 double qc_meda, qc_medb, qc_medc, qc_medq ;
481
482 /* Compute the QC params */
483 qc_meda = cpl_image_get_median(cpl_imagelist_get(fitres, 0)) ;
484 qc_medb = cpl_image_get_median(cpl_imagelist_get(fitres, 1)) ;
485 qc_medc = cpl_image_get_median(cpl_imagelist_get(fitres, 2)) ;
486 qc_medq = cpl_image_get_median(cpl_imagelist_get(fitres, 3)) ;
487
488 /* Get the QC params in qclist */
489 qclist = cpl_propertylist_new() ;
490 cpl_propertylist_append_double(qclist, "ESO QC DETLIN MEDA", qc_meda) ;
491 cpl_propertylist_append_double(qclist, "ESO QC DETLIN MEDB", qc_medb) ;
492 cpl_propertylist_append_double(qclist, "ESO QC DETLIN MEDC", qc_medc) ;
493 cpl_propertylist_append_double(qclist, "ESO QC DETLIN MEDQ", qc_medq) ;
494 cpl_propertylist_append_double(qclist, "ESO QC DETLIN LAMP",
495 sofi_img_detlin_config.lamp_stability ) ;
496
497 /* Write the _A FITS FILE */
498 irplib_dfs_save_image(set,
499 parlist,
500 set,
501 cpl_imagelist_get(fitres, 0),
502 CPL_BPP_IEEE_FLOAT,
503 "sofi_img_detlin",
504 SOFI_IMG_DETLIN_A,
505 qclist,
506 NULL,
507 PACKAGE "/" PACKAGE_VERSION,
508 "sofi_img_detlin_A.fits") ;
509
510 /* Write the _B FITS FILE */
511 irplib_dfs_save_image(set,
512 parlist,
513 set,
514 cpl_imagelist_get(fitres, 1),
515 CPL_BPP_IEEE_FLOAT,
516 "sofi_img_detlin",
517 SOFI_IMG_DETLIN_B,
518 qclist,
519 NULL,
520 PACKAGE "/" PACKAGE_VERSION,
521 "sofi_img_detlin_B.fits") ;
522
523 /* Write the _C FITS FILE */
524 irplib_dfs_save_image(set,
525 parlist,
526 set,
527 cpl_imagelist_get(fitres, 2),
528 CPL_BPP_IEEE_FLOAT,
529 "sofi_img_detlin",
530 SOFI_IMG_DETLIN_C,
531 qclist,
532 NULL,
533 PACKAGE "/" PACKAGE_VERSION,
534 "sofi_img_detlin_C.fits") ;
535
536 /* Write the _Q FITS FILE */
537 irplib_dfs_save_image(set,
538 parlist,
539 set,
540 cpl_imagelist_get(fitres, 3),
541 CPL_BPP_IEEE_FLOAT,
542 "sofi_img_detlin",
543 SOFI_IMG_DETLIN_Q,
544 qclist,
545 NULL,
546 PACKAGE "/" PACKAGE_VERSION,
547 "sofi_img_detlin_Q.fits") ;
548
549 /* Get the reference frame */
550 ref_frame = irplib_frameset_get_first_from_group(set, CPL_FRAME_GROUP_RAW) ;
551
552 /* Get FITS header from reference file */
553 if ((plist=cpl_propertylist_load(cpl_frame_get_filename(ref_frame),
554 0)) == NULL) {
555 cpl_msg_error(cpl_func, "getting header from reference frame");
556 cpl_propertylist_delete(qclist) ;
557 return CPL_ERROR_UNSPECIFIED;
558 }
559
560 /* Get the keywords for the paf file */
561 paflist = cpl_propertylist_new() ;
562 cpl_propertylist_copy_property_regexp(paflist, plist,
563 "^(ARCFILE|MJD-OBS|ESO TPL ID|DATE-OBS|ESO DET DIT|"
564 "ESO DET NDIT|ESO DET NCORRS|ESO DET MODE NAME)$", 0) ;
565 cpl_propertylist_delete(plist) ;
566
567 /* Copy the QC in paflist */
568 cpl_propertylist_copy_property_regexp(paflist, qclist, ".", 0) ;
569 cpl_propertylist_delete(qclist) ;
570
571 /* Save the PAF file */
572 cpl_dfs_save_paf("SOFI",
573 "sofi_img_detlin",
574 paflist,
575 "sofi_img_detlin_QC.paf") ;
576 cpl_propertylist_delete(paflist) ;
577 return 0;
578}
579
int sofi_dfs_set_groups(cpl_frameset *set)
Set the group as RAW or CALIB in a frameset.
Definition sofi_dfs.c:60
double sofi_pfits_get_dit(const cpl_propertylist *plist)
find out the DIT value
Definition sofi_pfits.c:195
const char * sofi_get_license(void)
Get the pipeline copyright and license.
Definition sofi_utils.c:60
cpl_frameset * sofi_extract_frameset(const cpl_frameset *in, const char *tag)
Extract the frames with the given tag from a frameset.
Definition sofi_utils.c:298