IIINSTRUMENT Pipeline Reference Manual 1.5.16
sofi_spc_jitter.c
1/* $Id: sofi_spc_jitter.c,v 1.45 2013-03-12 08:04:50 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:50 $
24 * $Revision: 1.45 $
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 <string.h>
37#include <math.h>
38#include <cpl.h>
39
40#include "irplib_plugin.h"
41#include "irplib_utils.h"
42#include "irplib_spectrum.h"
43
44#include "sofi_utils.h"
45#include "sofi_wavelength.h"
46#include "sofi_pfits.h"
47#include "sofi_dfs.h"
48
49/*-----------------------------------------------------------------------------
50 Defines
51 -----------------------------------------------------------------------------*/
52
53#define RECIPE_STRING "sofi_spc_jitter"
54
55#define SOFI_SPC_JITTER_OFFSET_ERR 10
56
57/*-----------------------------------------------------------------------------
58 Functions prototypes
59 -----------------------------------------------------------------------------*/
60
61static cpl_image ** sofi_spc_jitter_combine(cpl_frameset *, const char *,
62 const char *, const char *);
63static cpl_vector * sofi_spc_jitter_get_offsets(cpl_frameset *);
64static int * sofi_spc_jitter_classif(cpl_vector *, int *);
65static int off_comp(double, double, double);
66static cpl_imagelist * sofi_spc_jitter_saa_groups(cpl_imagelist *,
67 cpl_vector *, int *, int, cpl_vector **);
68static int sofi_spc_jitter_wavecal(const char *, const char *, cpl_image *,
69 cpl_frameset *);
70static cpl_imagelist * sofi_spc_jitter_nodded(cpl_imagelist *, cpl_vector *,
71 cpl_vector **);
72static cpl_imagelist * sofi_spc_jitter_distor(cpl_imagelist *, const char *);
73static double sofi_spc_jitter_refine_offset(cpl_image *, cpl_image *);
74static cpl_table * sofi_spc_jitter_extract(cpl_image *);
75static cpl_error_code sofi_spc_jitter_save(cpl_frameset *, const cpl_image *,
76 const cpl_table *,
77 const cpl_parameterlist *);
78
79cpl_recipe_define(sofi_spc_jitter, SOFI_BINARY_VERSION,
80 "Lars Lundin", PACKAGE_BUGREPORT, "2002,2003,2009",
81 "SOFI Spectro jitter recipe",
82 RECIPE_STRING " -- SOFI Spectro jitter recipe"
83 "The files listed in the Set Of Frames (sof-file) "
84 "must be tagged:\n"
85 "raw-file.fits "SOFI_SPC_JITTER_NODOBJ_RAW" or\n"
86 "raw-file.fits "SOFI_SPC_JITTER_NODSKY_RAW"\n"
87 "Calibration files:\n"
88 "oh-cat.fits "SOFI_CALPRO_OH_CAT" or\n"
89 "flat-file.fits "SOFI_CALIB_SPFLAT" or\n"
90 "arc-file.fits "SOFI_CALIB_ARC"\n");
91
92/*-----------------------------------------------------------------------------
93 Static variables
94 -----------------------------------------------------------------------------*/
95
96static struct {
97 /* Inputs */
98 int display;
99 int crosstalk;
100 /* wavecal_in : 0 for physical model, 1 for sky, 2 for arc table */
101 int wavecal_in;
102 int wavecal_degree;
103 double wavecal_err;
104 int wavecal_nsamples;
105 int wavecal_ppm;
106 int saa_refine;
107 double saa_rej_high;
108 double saa_rej_low;
109 int extr_spec_pos;
110 int extr_spec_width;
111 int extr_sky_le_width;
112 int extr_sky_ri_width;
113 int extr_sky_le_dist;
114 int extr_sky_ri_dist;
115 /* Outputs */
116 /* wavecal_out : 0 for physical model, 1 for sky, 2 for arc table */
117 cpl_vector * throws;
118 cpl_vector * intensities;
119 char filter[512];
120 int wavecal_out;
121 double wavecal_cc;
122 double wavecal_a0;
123 double wavecal_a1;
124 double wavecal_a2;
125 double wavecal_a3;
126 double wavecal_a4;
127} sofi_spc_jitter_config;
128
129
130/*-----------------------------------------------------------------------------
131 Functions code
132 -----------------------------------------------------------------------------*/
133
134/*----------------------------------------------------------------------------*/
142/*----------------------------------------------------------------------------*/
143static
144cpl_error_code sofi_spc_jitter_fill_parameterlist(cpl_parameterlist * self)
145{
146 const char * context = PACKAGE "." RECIPE_STRING;
147 cpl_error_code err;
148
149 cpl_ensure_code(self, CPL_ERROR_NULL_INPUT);
150
151 /* Fill the parameters list */
152
153
154 /* --crosstalk */
155 err = irplib_parameterlist_set_bool(self, PACKAGE, RECIPE_STRING,
156 "crosstalk", CPL_TRUE, NULL, context,
157 "Enable removal of the crosstalk "
158 "effect");
159 cpl_ensure_code(!err, err);
160
161 /* --wavecal */
162 err = irplib_parameterlist_set_string(self, PACKAGE, RECIPE_STRING,
163 "wavecal", "sky", NULL, context,
164 "Wavelength method: "
165 "phy or sky or arc");
166 cpl_ensure_code(!err, err);
167
168 /* --wavecal_ppm */
169 err = irplib_parameterlist_set_bool(self, PACKAGE, RECIPE_STRING,
170 "wavecal_ppm", CPL_FALSE, NULL, context,
171 "Enable Point Pattern Matching in the "
172 "wavelength calibration");
173 cpl_ensure_code(!err, err);
174
175 /* --wavecal_degree */
176 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
177 "wavecal_degree", 2, NULL, context,
178 "Degree of the wavelength dispersion "
179 "polynomial");
180 cpl_ensure_code(!err, err);
181
182 /* --wavecal_nsamples */
183 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
184 "wavecal_nsamples", 100, NULL, context,
185 "Number of samples for the wavelength "
186 "calibration");
187 cpl_ensure_code(!err, err);
188
189 /* --wavecal_err */
190 err = irplib_parameterlist_set_double(self, PACKAGE, RECIPE_STRING,
191 "wavecal_err", 1000.0, NULL, context,
192 "The wavelength error [Angstrom]");
193 cpl_ensure_code(!err, err);
194
195 /* --saa_refine */
196 err = irplib_parameterlist_set_bool(self, PACKAGE, RECIPE_STRING,
197 "saa_refine", CPL_TRUE, NULL, context,
198 "Enable refinement of the offsets");
199 cpl_ensure_code(!err, err);
200
201 /* --saa_rej */
202 err = irplib_parameterlist_set_string(self, PACKAGE, RECIPE_STRING,
203 "saa_rej", "0.1,0.1", NULL, context,
204 "Low, high SAA rejection [%,%]");
205 cpl_ensure_code(!err, err);
206
207 /* --spec_pos */
208 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
209 "spec_pos", -1, NULL, context,
210 "Spectrum position");
211 cpl_ensure_code(!err, err);
212
213 /* --spec_width */
214 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
215 "spec_width", 10, NULL, context,
216 "Spectrum width");
217 cpl_ensure_code(!err, err);
218
219 /* --sky_le_width */
220 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
221 "sky_le_width", 10, NULL, context,
222 "Sky width left of the spectrum");
223 cpl_ensure_code(!err, err);
224
225 /* --sky_ri_width */
226 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
227 "sky_ri_width", 10, NULL, context,
228 "Sky width right of the spectrum");
229 cpl_ensure_code(!err, err);
230
231 /* --sky_le_dist */
232 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
233 "sky_le_dist", -1, NULL, context,
234 "Sky distance left of the spectrum");
235 cpl_ensure_code(!err, err);
236
237 /* --sky_ri_dist */
238 err = irplib_parameterlist_set_int(self, PACKAGE, RECIPE_STRING,
239 "sky_ri_dist", -1, NULL, context,
240 "Sky distance right of the spectrum");
241 cpl_ensure_code(!err, err);
242
243 /* --display */
244 err = irplib_parameterlist_set_bool(self, PACKAGE, RECIPE_STRING, "display",
245 CPL_FALSE, NULL, context,
246 "Enable plotting");
247 cpl_ensure_code(!err, err);
248
249 return CPL_ERROR_NONE;
250}
251
252/*----------------------------------------------------------------------------*/
259/*----------------------------------------------------------------------------*/
260static int sofi_spc_jitter(cpl_frameset * framelist,
261 const cpl_parameterlist * parlist)
262{
263 cpl_propertylist * plist;
264 const char * sval;
265 cpl_frameset * rawframes = NULL;
266 const cpl_frame * cur_frame;
267 const char * flat;
268 const char * oh;
269 const char * arc;
270 cpl_image ** combined = NULL;
271 cpl_table * extracted = NULL;
272
273 /* Initialise */
274 sofi_spc_jitter_config.wavecal_out = -1;
275 sofi_spc_jitter_config.wavecal_cc = -1.0;
276 sofi_spc_jitter_config.throws = NULL;
277 sofi_spc_jitter_config.intensities = NULL;
278 sofi_spc_jitter_config.filter[0] = (char)0;
279
280 /* Retrieve input parameters */
281
282 /* Cross-talk */
283 sofi_spc_jitter_config.crosstalk
284 = irplib_parameterlist_get_bool(parlist, PACKAGE, RECIPE_STRING,
285 "crosstalk");
286
287 /* --wavecal */
288 sval = irplib_parameterlist_get_string(parlist, PACKAGE, RECIPE_STRING,
289 "wavecal");
290 bug_if(sval == NULL);
291 if (!strcmp(sval, "phy")) sofi_spc_jitter_config.wavecal_in = 0;
292 else if (!strcmp(sval, "sky")) sofi_spc_jitter_config.wavecal_in = 1;
293 else if (!strcmp(sval, "arc")) sofi_spc_jitter_config.wavecal_in = 2;
294 else {
295 error_if(1, CPL_ERROR_UNSUPPORTED_MODE,
296 "Invalid value for wavecal option: %s", sval);
297 }
298
299 /* --wavecal_degree */
300 sofi_spc_jitter_config.wavecal_degree
301 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
302 "wavecal_degree");
303
304 /* --wavecal_nsamples */
305 sofi_spc_jitter_config.wavecal_nsamples
306 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
307 "wavecal_nsamples");
308
309 /* --wavecal_err */
310 sofi_spc_jitter_config.wavecal_err
311 = irplib_parameterlist_get_double(parlist, PACKAGE, RECIPE_STRING,
312 "wavecal_err");
313
314 /* --wavecal_ppm */
315 sofi_spc_jitter_config.wavecal_ppm
316 = irplib_parameterlist_get_bool(parlist, PACKAGE, RECIPE_STRING,
317 "wavecal_ppm");
318
319 /* Refine of offsets */
320 sofi_spc_jitter_config.saa_refine
321 = irplib_parameterlist_get_bool(parlist, PACKAGE, RECIPE_STRING,
322 "saa_refine");
323
324 /* Rejection parameters for SAA */
325 sval = irplib_parameterlist_get_string(parlist, PACKAGE, RECIPE_STRING,
326 "saa_rej");
327 skip_if (sscanf(sval, "%lg,%lg",
328 &sofi_spc_jitter_config.saa_rej_low,
329 &sofi_spc_jitter_config.saa_rej_high) != 2);
330
331 /* --spec_pos */
332 sofi_spc_jitter_config.extr_spec_pos
333 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
334 "spec_pos");
335
336 /* --spec_width */
337 sofi_spc_jitter_config.extr_spec_width
338 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
339 "spec_width");
340
341 /* --sky_le_width */
342 sofi_spc_jitter_config.extr_sky_le_width
343 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
344 "sky_le_width");
345 /* --sky_ri_width */
346 sofi_spc_jitter_config.extr_sky_ri_width
347 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
348 "sky_ri_width");
349
350 /* --sky_le_dist */
351 sofi_spc_jitter_config.extr_sky_le_dist
352 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
353 "sky_le_dist");
354
355 /* --sky_ri_dist */
356 sofi_spc_jitter_config.extr_sky_ri_dist
357 = irplib_parameterlist_get_int(parlist, PACKAGE, RECIPE_STRING,
358 "sky_ri_dist");
359 /* Display */
360 sofi_spc_jitter_config.display
361 = irplib_parameterlist_get_bool(parlist, PACKAGE, RECIPE_STRING,
362 "display");
363
364 /* Identify the RAW and CALIB frames in the input frameset */
365 skip_if (sofi_dfs_set_groups(framelist));
366
367 /* Get the filter */
368 cur_frame = cpl_frameset_get_position(framelist, 0);
369 plist = cpl_propertylist_load(cpl_frame_get_filename(cur_frame), 0);
370 skip_if(plist == NULL);
371
372 if ((sval = sofi_pfits_get_filter(plist)) != NULL) {
373 strcpy(sofi_spc_jitter_config.filter, sval); /* FIXME: Replace */
374 } else {
375 cpl_error_reset();
376 }
377 cpl_propertylist_delete(plist);
378
379 /* Retrieve calibration data */
380 arc = sofi_extract_filename(framelist, SOFI_CALIB_ARC);
381 flat = sofi_extract_filename(framelist, SOFI_CALIB_SPFLAT);
382 oh = sofi_extract_filename(framelist, SOFI_CALPRO_OH_CAT);
383
384 /* Retrieve raw frames */
385 rawframes = sofi_extract_frameset(framelist, SOFI_SPC_JITTER_NODOBJ_RAW);
386 error_if (rawframes == NULL, CPL_ERROR_DATA_NOT_FOUND,
387 "Cannot find the raw frames in the input list");
388
389 /* If arc calibration requested, arc file must be provided */
390 error_if (arc == NULL && sofi_spc_jitter_config.wavecal_in == 2,
391 CPL_ERROR_ILLEGAL_INPUT,
392 "You must provide an ARC file for the wl cal");
393
394 /* Create the combined image */
395 cpl_msg_info(cpl_func, "Create the combined image");
396 cpl_msg_indent_more();
397 combined = sofi_spc_jitter_combine(rawframes, oh, flat, arc);
398 cpl_msg_indent_less();
399
400 error_if (combined==NULL, CPL_ERROR_ILLEGAL_INPUT,
401 "Cannot combine the images");
402
403 /* Extract the spectrum */
404 cpl_msg_info(cpl_func, "Extract the spectrum");
405 cpl_msg_indent_more();
406 extracted = sofi_spc_jitter_extract(combined[0]);
407 cpl_msg_indent_less();
408
409 error_if (extracted == NULL, CPL_ERROR_DATA_NOT_FOUND,
410 "Cannot extract the spectrum");
411
412 /* Write the products */
413 cpl_msg_info(cpl_func, "Save the products");
414 skip_if(sofi_spc_jitter_save(framelist, combined[0], extracted, parlist));
415
416 end_skip;
417
418 if (combined != NULL) {
419 cpl_image_delete(combined[0]);
420 cpl_image_delete(combined[1]);
421 cpl_free(combined);
422 }
423
424 cpl_vector_delete(sofi_spc_jitter_config.intensities);
425 cpl_vector_delete(sofi_spc_jitter_config.throws);
426
427 cpl_frameset_delete(rawframes);
428
429 cpl_table_delete(extracted);
430
431 return cpl_error_get_code();
432}
433
434/*----------------------------------------------------------------------------*/
443/*----------------------------------------------------------------------------*/
444static cpl_image ** sofi_spc_jitter_combine(
445 cpl_frameset * rawframes,
446 const char * oh,
447 const char * flat,
448 const char * arc)
449{
450 cpl_imagelist * ilist;
451 cpl_image * tmp_im;
452 cpl_vector * offsets;
453 int * groups;
454 int ngroups;
455 cpl_imagelist * abba;
456 cpl_vector * abba_off;
457 cpl_imagelist * nodded;
458 cpl_vector * nodded_off_x;
459 cpl_vector * nodded_off_y;
460 double throw;
461 cpl_table * extracted;
462 double intensity;
463 double * pnodded_off_x;
464 cpl_imagelist * nodded_warped;
465 cpl_bivector * nodded_offsets;
466 cpl_image ** combined;
467 int nima;
468 double new_offset;
469 int i;
470
471 /* Test entries */
472 if (rawframes == NULL) return NULL;
473
474 /* Load the images */
475 cpl_msg_info(cpl_func, "Load the data");
476 cpl_msg_indent_more();
477 if ((ilist = cpl_imagelist_load_frameset(rawframes, CPL_TYPE_FLOAT, 1,
478 0)) == NULL) {
479 cpl_msg_error(cpl_func, "cannot load the data");
480 cpl_msg_indent_less();
481 return NULL;
482 }
483 cpl_msg_indent_less();
484
485 /* Apply the cross-talk correction */
486 if (sofi_spc_jitter_config.crosstalk) {
487 cpl_msg_info(cpl_func, "Apply the cross-talk correction");
488 cpl_msg_indent_more();
489 if (sofi_correct_crosstalk_list(ilist) == -1) {
490 cpl_msg_error(cpl_func, "Cannot correct for Cross-talk");
491 cpl_imagelist_delete(ilist);
492 cpl_msg_indent_less();
493 return NULL;
494 }
495 cpl_msg_indent_less();
496 }
497
498 /* Apply the flafield */
499 if (flat != NULL) {
500 cpl_msg_info(cpl_func, "Apply the flatfield correction");
501 if ((tmp_im = cpl_image_load(flat, CPL_TYPE_FLOAT, 0, 0)) == NULL) {
502 cpl_msg_warning(cpl_func, "cannot load the flat field");
503 } else {
504 if (cpl_imagelist_divide_image(ilist, tmp_im) != CPL_ERROR_NONE) {
505 cpl_msg_warning(cpl_func, "cannot apply the flat field");
506 }
507 cpl_image_delete(tmp_im);
508 }
509 }
510
511 /* Get the offsets */
512 cpl_msg_info(cpl_func, "Get the offsets");
513 if ((offsets = sofi_spc_jitter_get_offsets(rawframes)) == NULL) {
514 cpl_msg_error(cpl_func, "cannot get the offsets");
515 cpl_imagelist_delete(ilist);
516 return NULL;
517 }
518
519 /* Classify in groups the a and b images sequence */
520 cpl_msg_info(cpl_func, "Classify in groups");
521 cpl_msg_indent_more();
522 if ((groups = sofi_spc_jitter_classif(offsets, &ngroups)) == NULL) {
523 cpl_msg_error(cpl_func, "cannot classify the data");
524 cpl_imagelist_delete(ilist);
525 cpl_vector_delete(offsets);
526 cpl_msg_indent_less();
527 return NULL;
528 }
529 cpl_msg_indent_less();
530
531 /* Shift and add each group to one image */
532 cpl_msg_info(cpl_func, "Shift and add each group to one image");
533 cpl_msg_indent_more();
534 if ((abba = sofi_spc_jitter_saa_groups(ilist, offsets, groups,
535 ngroups, &abba_off)) == NULL) {
536 cpl_msg_error(cpl_func, "cannot shift and add groups");
537 cpl_imagelist_delete(ilist);
538 cpl_vector_delete(offsets);
539 cpl_free(groups);
540 cpl_msg_indent_less();
541 return NULL;
542 }
543 cpl_imagelist_delete(ilist);
544 cpl_free(groups);
545 cpl_vector_delete(offsets);
546 cpl_msg_indent_less();
547
548 /* Compute the wavelength calibration */
549 cpl_msg_info(cpl_func, "Compute the wavelength calibration");
550 cpl_msg_indent_more();
551 if (sofi_spc_jitter_wavecal(arc, oh, cpl_imagelist_get(abba, 0),
552 rawframes) == -1) {
553 cpl_msg_error(cpl_func, "cannot compute the wavelength");
554 cpl_imagelist_delete(abba);
555 cpl_vector_delete(abba_off);
556 cpl_msg_indent_less();
557 return NULL;
558 }
559 cpl_msg_indent_less();
560
561 /* Create the nodded images */
562 cpl_msg_info(cpl_func, "Create the nodded images");
563 cpl_msg_indent_more();
564 if ((nodded = sofi_spc_jitter_nodded(abba, abba_off,
565 &nodded_off_x))==NULL) {
566 cpl_msg_error(cpl_func, "cannot create the nodded images");
567 cpl_imagelist_delete(abba);
568 cpl_vector_delete(abba_off);
569 cpl_msg_indent_less();
570 return NULL;
571 }
572 cpl_imagelist_delete(abba);
573 cpl_msg_indent_less();
574
575 /* Get the throw offsets from abba_off */
576 nima = cpl_imagelist_get_size(nodded);
577 sofi_spc_jitter_config.throws = cpl_vector_new(nima);
578 for (i=0; i<nima/2; i++) {
579 throw = fabs( (cpl_vector_get(abba_off, 2*i))-
580 (cpl_vector_get(abba_off, 2*i+1)));
581 cpl_vector_set(sofi_spc_jitter_config.throws, 2*i, throw);
582 cpl_vector_set(sofi_spc_jitter_config.throws, 2*i+1, throw);
583 }
584 cpl_vector_delete(abba_off);
585
586 /* Get the spectra intensities in the nodded images */
587 cpl_msg_info(cpl_func, "Compute the spectra intensities");
588 cpl_msg_indent_more();
589 nima = cpl_imagelist_get_size(nodded);
590 sofi_spc_jitter_config.intensities = cpl_vector_new(nima);
591 for (i=0; i<nima; i++) {
592 if ((extracted = sofi_spc_jitter_extract(
593 cpl_imagelist_get(nodded, i))) == NULL) {
594 cpl_msg_warning(cpl_func,"Cannot extract the spectrum from nodded %d",
595 i+1);
596 intensity = -1.0;
597 } else {
598 intensity = cpl_table_get_column_mean(extracted,
599 "Extracted_spectrum_value");
600 intensity *= cpl_table_get_nrow(extracted);
601 cpl_table_delete(extracted);
602 }
603 cpl_msg_info(cpl_func, "Spectrum intensity nb %d: %g", i+1,
604 intensity);
605 cpl_vector_set(sofi_spc_jitter_config.intensities, i, intensity);
606 }
607 cpl_msg_indent_less();
608
609 /* Distortion correction */
610 if (arc) {
611 cpl_msg_info(cpl_func, "Correct the distortion on nodded images");
612 cpl_msg_indent_more();
613 if ((nodded_warped = sofi_spc_jitter_distor(nodded, arc)) == NULL) {
614 cpl_msg_error(cpl_func, "cannot correct the distortion");
615 cpl_imagelist_delete(nodded);
616 cpl_vector_delete(nodded_off_x);
617 cpl_msg_indent_less();
618 return NULL;
619 }
620 cpl_imagelist_delete(nodded);
621 nodded = nodded_warped;
622 cpl_msg_indent_less();
623 }
624
625 /* Refine the offsets if requested */
626 if (sofi_spc_jitter_config.saa_refine) {
627 cpl_msg_info(cpl_func, "Refine the offsets");
628 pnodded_off_x = cpl_vector_get_data(nodded_off_x);
629 for (i=0; i<cpl_imagelist_get_size(nodded); i++) {
630 new_offset = sofi_spc_jitter_refine_offset(
631 cpl_imagelist_get(nodded, 0),
632 cpl_imagelist_get(nodded, i));
633 if (new_offset > 5000) {
634 cpl_msg_debug(cpl_func, "cannot refine the offset - keep %g",
635 pnodded_off_x[i]);
636 } else {
637 if (fabs(new_offset-pnodded_off_x[i]) <
638 SOFI_SPC_JITTER_OFFSET_ERR) {
639 cpl_msg_debug(cpl_func, "refined offset : %g (old was %g)",
640 new_offset, pnodded_off_x[i]);
641 pnodded_off_x[i] = new_offset;
642 } else {
643 cpl_msg_debug(cpl_func,
644 "refined offset %g too different - keep %g",
645 new_offset, pnodded_off_x[i]);
646 }
647 }
648 }
649 }
650
651 /* Images combination */
652 /* Get the offsets in a bivector */
653 nodded_off_y = cpl_vector_duplicate(nodded_off_x);
654 cpl_vector_fill(nodded_off_y, 0.0);
655 nodded_offsets = cpl_bivector_wrap_vectors(nodded_off_x, nodded_off_y);
656 /* Shift and add */
657 cpl_msg_info(cpl_func, "Apply the shift and add on the nodded frames");
658 nima = cpl_imagelist_get_size(nodded);
659 if ((combined = cpl_geom_img_offset_saa(nodded, nodded_offsets,
660 CPL_KERNEL_DEFAULT,
661 (int)(sofi_spc_jitter_config.saa_rej_low * nima),
662 (int)(sofi_spc_jitter_config.saa_rej_high * nima),
663 CPL_GEOM_FIRST, NULL, NULL)) == NULL) {
664 cpl_msg_error(cpl_func, "Cannot shift and add group");
665 cpl_imagelist_delete(nodded);
666 cpl_bivector_unwrap_vectors(nodded_offsets);
667 cpl_vector_delete(nodded_off_x);
668 cpl_vector_delete(nodded_off_y);
669 return NULL;
670 }
671 cpl_imagelist_delete(nodded);
672 cpl_bivector_unwrap_vectors(nodded_offsets);
673 cpl_vector_delete(nodded_off_x);
674 cpl_vector_delete(nodded_off_y);
675 return combined;
676}
677
678/*----------------------------------------------------------------------------*/
684/*----------------------------------------------------------------------------*/
685static cpl_vector * sofi_spc_jitter_get_offsets(cpl_frameset * rawframes)
686{
687 cpl_vector * offsets;
688 double * pvect;
689 int nraw;
690 cpl_frame * cur_frame;
691 cpl_propertylist * plist;
692 int i;
693
694 /* Test entries */
695 if (rawframes == NULL) return NULL;
696
697 /* Initialise */
698 nraw = cpl_frameset_get_size(rawframes);
699
700 /* Get the rawframes X offsets */
701 offsets = cpl_vector_new(nraw);
702 pvect = cpl_vector_get_data(offsets);
703 for (i=0; i<nraw; i++) {
704 cur_frame = cpl_frameset_get_position(rawframes, i);
705 if ((plist=cpl_propertylist_load(cpl_frame_get_filename(cur_frame),
706 0)) == NULL) {
707 cpl_msg_error(cpl_func, "cannot get property list");
708 cpl_vector_delete(offsets);
709 return NULL;
710 }
711 pvect[i] = sofi_pfits_get_cumoffsetx(plist);
712 if (cpl_error_get_code()) {
713 cpl_msg_error(cpl_func, "cannot get the offset from the header");
714 cpl_vector_delete(offsets);
715 cpl_propertylist_delete(plist);
716 return NULL;
717 }
718 cpl_propertylist_delete(plist);
719 pvect[i] *= -1.0;
720 }
721 return offsets;
722}
723
724/*----------------------------------------------------------------------------*/
762/*----------------------------------------------------------------------------*/
763static int * sofi_spc_jitter_classif(
764 cpl_vector * offsets,
765 int * ngroups)
766{
767 double * pvect;
768 int nraw;
769 double offset_thresh;
770 cpl_vector * tmp_vec;
771 int * groups;
772 int last_group;
773 int i, j, k, l;
774
775 /* Test entries */
776 if (offsets == NULL) return NULL;
777
778 /* Initialise */
779 nraw = cpl_vector_get_size(offsets);
780
781 /* Separate the offsets in 2 categories */
782 tmp_vec = cpl_vector_duplicate(offsets);
783 cpl_vector_sort(tmp_vec, 1);
784 pvect = cpl_vector_get_data(tmp_vec);
785 if (pvect[0] == pvect[nraw-1]) {
786 cpl_msg_error(cpl_func, "Only one offset in the list - abort");
787 cpl_vector_delete(tmp_vec);
788 return NULL;
789 }
790 offset_thresh = (pvect[0] + pvect[nraw-1]) / 2.0;
791 cpl_vector_delete(tmp_vec);
792
793 /* Identify the different A and B groups */
794 pvect = cpl_vector_get_data(offsets);
795 *ngroups = 0;
796 groups = cpl_calloc(nraw, sizeof(int));
797
798 /* Create a look up table to associate the ith obj with the jth frame */
799 i = 0;
800 while (i < nraw) {
801 j = 0;
802 /* Count the number of successive '+' or '-' (j) */
803 while ((i+j<nraw) &&
804 (!off_comp(pvect[i], pvect[i+j], offset_thresh))) j++;
805
806 if (i+j >= nraw) i = nraw;
807 else {
808 k = 0;
809 /* Check if there are j '-' or '+' (k) */
810 while ((i+j+k < nraw)
811 && (!off_comp(pvect[i+j], pvect[i+j+k], offset_thresh))
812 && (k<j)) k++;
813 last_group = 1;
814 if (i+j+k < nraw) {
815 for (l=i+j+k; l<nraw; l++) {
816 if (off_comp(pvect[i+j], pvect[l], offset_thresh)) {
817 last_group = 0;
818 break;
819 }
820 }
821 }
822 if (last_group == 0) {
823 for (l=0; l<j; l++) groups[i+l] = *ngroups + 1;
824 for (l=0; l<k; l++) groups[i+j+l] = *ngroups + 2;
825 *ngroups += 2;
826 i += j+k;
827 } else {
828 for (l=0; l<j; l++) groups[i+l] = *ngroups + 1;
829 for (l=0; l<nraw - (i+j); l++) groups[i+j+l] =*ngroups + 2;
830 *ngroups += 2;
831 i = nraw;
832 }
833 }
834 }
835
836 /* Nb of groups found should be even */
837 if (*ngroups % 2) {
838 cpl_msg_error(cpl_func, "Odd number of groups found");
839 cpl_free(groups);
840 return NULL;
841 }
842
843 return groups;
844}
845
846/*----------------------------------------------------------------------------*/
877/*----------------------------------------------------------------------------*/
878static cpl_imagelist * sofi_spc_jitter_saa_groups(
879 cpl_imagelist * ilist,
880 cpl_vector * offsets,
881 int * groups,
882 int ngroups,
883 cpl_vector ** abba_off)
884{
885 cpl_imagelist * abba;
886 cpl_imagelist * group_list;
887 cpl_image * tmp_ima;
888 cpl_image ** combined;
889 cpl_bivector * group_off;
890 double * pgroup_off;
891 double * poffsets;
892 double * pabba_off;
893 int nima;
894 int saa;
895 int i, j, k;
896
897 /* Test entries */
898 if ((ilist == NULL) || (offsets == NULL) || (groups == NULL)) return NULL;
899
900 /* Initialise */
901 nima = cpl_imagelist_get_size(ilist);
902 poffsets = cpl_vector_get_data(offsets);
903
904 /* Create the output image list */
905 abba = cpl_imagelist_new();
906 *abba_off = cpl_vector_new(ngroups);
907 pabba_off = cpl_vector_get_data(*abba_off);
908
909 /* Loop on the groups */
910 for (i=0; i<ngroups; i++) {
911 /* Initialise */
912 saa = 0;
913 /* Create the group list of images */
914 group_list = cpl_imagelist_new();
915 k = 0;
916 for (j=0; j<nima; j++) {
917 if (i+1 == groups[j]) {
918 /* Get the first offset of the group in abba_off */
919 if (k==0) pabba_off[i] = poffsets[j];
920 /* To know if we need the saa (shift and add) */
921 if (fabs(pabba_off[i]-poffsets[j]) > 1e-3) saa = 1;
922 /* Copy the images of the group in group_list */
923 tmp_ima = cpl_image_duplicate(cpl_imagelist_get(ilist, j));
924 cpl_imagelist_set(group_list, tmp_ima, k);
925 tmp_ima = NULL;
926 k++;
927 }
928 }
929
930 if (saa) {
931 /* Get the offsets of the group in group_off */
932 group_off = cpl_bivector_new(k);
933 cpl_vector_fill(cpl_bivector_get_y(group_off), 0.0);
934 pgroup_off = cpl_bivector_get_x_data(group_off);
935 k = 0;
936 for (j=0; j<nima; j++) {
937 if (i+1 == groups[j]) {
938 pgroup_off[k] = poffsets[j];
939 k++;
940 }
941 }
942 cpl_vector_subtract_scalar(cpl_bivector_get_x(group_off),
943 pabba_off[i]);
944 /* Shift and add */
945 cpl_msg_debug(cpl_func, "Apply shift-and-add for group %d", i+1);
946 if ((combined = cpl_geom_img_offset_saa(group_list,
947 group_off, CPL_KERNEL_DEFAULT, 0, 0,
948 CPL_GEOM_FIRST, NULL, NULL)) == NULL) {
949 cpl_msg_error(cpl_func,"Cannot shift and add group nb %d", i+1);
950 cpl_imagelist_delete(group_list);
951 cpl_bivector_delete(group_off);
952 cpl_imagelist_delete(abba);
953 cpl_vector_delete(*abba_off);
954 return NULL;
955 }
956 cpl_bivector_delete(group_off);
957 cpl_image_delete(combined[1]);
958 cpl_imagelist_set(abba, combined[0], i);
959 cpl_free(combined);
960 } else {
961 /* Averaging */
962 cpl_msg_debug(cpl_func, "Apply averaging for group %d", i+1);
963 if ((tmp_ima = cpl_imagelist_collapse_create(group_list)) == NULL) {
964 cpl_msg_error(cpl_func, "Cannot average group nb %d", i+1);
965 cpl_imagelist_delete(group_list);
966 cpl_imagelist_delete(abba);
967 cpl_vector_delete(*abba_off);
968 return NULL;
969 }
970 cpl_imagelist_set(abba, tmp_ima, i);
971 }
972 cpl_imagelist_delete(group_list);
973 }
974 return abba;
975}
976
977/*----------------------------------------------------------------------------*/
986/*----------------------------------------------------------------------------*/
987static int sofi_spc_jitter_wavecal(
988 const char * arc,
989 const char * oh,
990 cpl_image * ima,
991 cpl_frameset * raw)
992{
993 cpl_table * arc_tab;
994 cpl_polynomial * phdisprel;
995 cpl_frame * cur_frame;
996 const char * cur_fname;
997 cpl_polynomial * disprel;
998 double slit_width;
999 cpl_size p[5];
1000 double xc, a, b, c, d, e;
1001 int degree;
1002
1003 /* Get the wavelength from the arc file */
1004 if (sofi_spc_jitter_config.wavecal_in == 2) {
1005 if (arc == NULL) {
1006 cpl_msg_error(cpl_func, "Missing arc for the wavelength calib");
1007 return CPL_ERROR_UNSPECIFIED;
1008 }
1009 cpl_msg_info(cpl_func, "Get the wavelength from the ARC file");
1010 if ((arc_tab = cpl_table_load(arc, 1, 0)) == NULL) {
1011 cpl_msg_error(cpl_func, "Cannot load the arc table");
1012 sofi_spc_jitter_config.wavecal_out = -1;
1013 return CPL_ERROR_UNSPECIFIED;
1014 }
1015 sofi_spc_jitter_config.wavecal_a0 =
1016 cpl_table_get_double(arc_tab, "WL_coefficients", 0, NULL);
1017 sofi_spc_jitter_config.wavecal_a1 =
1018 cpl_table_get_double(arc_tab, "WL_coefficients", 1, NULL);
1019 sofi_spc_jitter_config.wavecal_a2 =
1020 cpl_table_get_double(arc_tab, "WL_coefficients", 2, NULL);
1021 sofi_spc_jitter_config.wavecal_a3 =
1022 cpl_table_get_double(arc_tab, "WL_coefficients", 3, NULL);
1023 sofi_spc_jitter_config.wavecal_a4 =
1024 cpl_table_get_double(arc_tab, "WL_coefficients", 4, NULL);
1025 cpl_table_delete(arc_tab);
1026 sofi_spc_jitter_config.wavecal_out = 2;
1027 sofi_spc_jitter_config.wavecal_cc = -1.0;
1028 return 0;
1029 }
1030
1031 /* Get the reference frame */
1032 cur_frame = cpl_frameset_get_position(raw, 0);
1033 cur_fname = cpl_frame_get_filename(cur_frame);
1034
1035 /* Get the physical model */
1036 cpl_msg_info(cpl_func, "Compute the physical model");
1037 cpl_msg_indent_more();
1038 if ((phdisprel = sofi_get_disprel_estimate(cur_fname, 1)) == NULL) {
1039 cpl_msg_error(cpl_func, "cannot compute the physical model");
1040 sofi_spc_jitter_config.wavecal_out = -1;
1041 cpl_msg_indent_less();
1042 return CPL_ERROR_UNSPECIFIED;
1043 }
1044 p[0] = 0; p[1] = 1;
1045 a = cpl_polynomial_get_coeff(phdisprel, p);
1046 b = cpl_polynomial_get_coeff(phdisprel, p + 1);
1047 cpl_msg_info(cpl_func, "f(x)=%g + %g*x", a, b);
1048 sofi_spc_jitter_config.wavecal_a0 = a;
1049 sofi_spc_jitter_config.wavecal_a1 = b;
1050 sofi_spc_jitter_config.wavecal_a2 = 0.0;
1051 sofi_spc_jitter_config.wavecal_a3 = 0.0;
1052 sofi_spc_jitter_config.wavecal_a4 = 0.0;
1053 sofi_spc_jitter_config.wavecal_cc = -1.0;
1054 sofi_spc_jitter_config.wavecal_out = 0;
1055 cpl_msg_indent_less();
1056
1057 /* Compute the wavelength using the sky lines */
1058 if (sofi_spc_jitter_config.wavecal_in == 1) {
1059 /* Compute the slit_width */
1060 if ((slit_width = sofi_get_slitwidth(cur_fname)) == -1) {
1061 cpl_msg_warning(cpl_func, "cannot get the slit width");
1062 cpl_polynomial_delete(phdisprel);
1063 return 0;
1064 }
1065 /* Compute the wavelength */
1066 cpl_msg_info(cpl_func, "Compute the wavelength with the sky lines");
1067 cpl_msg_indent_more();
1068 if ((disprel = sofi_wavelength_engine(ima, "oh", oh, NULL, NULL,
1069 phdisprel, slit_width,
1070 sofi_spc_jitter_config.wavecal_degree,
1071 sofi_spc_jitter_config.wavecal_err,
1072 sofi_spc_jitter_config.wavecal_nsamples,
1073 sofi_spc_jitter_config.wavecal_ppm,
1074 sofi_spc_jitter_config.display, &xc)) == NULL) {
1075 cpl_msg_error(cpl_func, "cannot compute the dispersion relation");
1076 cpl_polynomial_delete(phdisprel);
1077 cpl_msg_indent_less();
1078 return 0;
1079 }
1080 cpl_msg_info(cpl_func, "Cross correlation factor: %g", xc);
1081 p[0] = 0; p[1] = 1; p[2] = 2; p[3] = 3; p[4] = 4;
1082 a = b = c = d = e = 0.0;
1083 degree = cpl_polynomial_get_degree(disprel);
1084 a = cpl_polynomial_get_coeff(disprel, p);
1085 if (degree > 0) b = cpl_polynomial_get_coeff(disprel, p + 1);
1086 if (degree > 1) c = cpl_polynomial_get_coeff(disprel, p + 2);
1087 if (degree > 2) d = cpl_polynomial_get_coeff(disprel, p + 3);
1088 if (degree > 3) e = cpl_polynomial_get_coeff(disprel, p + 4);
1089 if (degree == 0)
1090 cpl_msg_info(cpl_func, "f(x)=%g", a);
1091 if (degree == 1)
1092 cpl_msg_info(cpl_func, "f(x)=%g + %g*x", a, b);
1093 if (degree == 2)
1094 cpl_msg_info(cpl_func, "f(x)=%g + %g*x + %g*x^2", a, b, c);
1095 if (degree == 3)
1096 cpl_msg_info(cpl_func, "f(x)=%g + %g*x + %g*x^2 + %g*x^3",
1097 a, b, c, d);
1098 if (degree == 4)
1099 cpl_msg_info(cpl_func, "f(x)=%g + %g*x + %g*x^2 + %g*x^3 + %g*x^4",
1100 a, b, c, d, e);
1101 sofi_spc_jitter_config.wavecal_a0 = a;
1102 sofi_spc_jitter_config.wavecal_a1 = b;
1103 sofi_spc_jitter_config.wavecal_a2 = c;
1104 sofi_spc_jitter_config.wavecal_a3 = d;
1105 sofi_spc_jitter_config.wavecal_a4 = e;
1106 sofi_spc_jitter_config.wavecal_cc = xc;
1107 sofi_spc_jitter_config.wavecal_out = 1;
1108 cpl_polynomial_delete(disprel);
1109 cpl_msg_indent_less();
1110 }
1111 cpl_polynomial_delete(phdisprel);
1112 return 0;
1113}
1114
1115/*----------------------------------------------------------------------------*/
1136/*----------------------------------------------------------------------------*/
1137static cpl_imagelist * sofi_spc_jitter_nodded(
1138 cpl_imagelist * abba,
1139 cpl_vector * abba_off,
1140 cpl_vector ** nodded_off)
1141{
1142 cpl_imagelist * nodded;
1143 cpl_image * tmp_ima;
1144 int nima;
1145 double ref_off;
1146 int i;
1147
1148 /* Test entries */
1149 if ((abba == NULL) || (abba_off == NULL)) return NULL;
1150
1151 /* Initialise */
1152 nima = cpl_imagelist_get_size(abba);
1153 if (nima % 2) {
1154 cpl_msg_error(cpl_func, "Number of images should be even");
1155 return NULL;
1156 }
1157
1158 /* Create the offsets between the nodded images */
1159 *nodded_off = cpl_vector_duplicate(abba_off);
1160 /* The image list to contain the nodded images */
1161 nodded = cpl_imagelist_new();
1162 for (i=0; i<(nima/2); i++) {
1163 /* a-b */
1164 tmp_ima = cpl_image_duplicate(cpl_imagelist_get(abba, 2*i));
1165 cpl_image_subtract(tmp_ima, cpl_imagelist_get(abba, 2*i+1));
1166 cpl_imagelist_set(nodded, tmp_ima, 2*i);
1167 /* b-a */
1168 tmp_ima = cpl_image_duplicate(cpl_imagelist_get(abba, 2*i+1));
1169 cpl_image_subtract(tmp_ima, cpl_imagelist_get(abba, 2*i));
1170 cpl_imagelist_set(nodded, tmp_ima, 2*i+1);
1171 }
1172
1173 /* Subtract the first offset to the others */
1174 ref_off = cpl_vector_get(*nodded_off, 0);
1175 cpl_vector_subtract_scalar(*nodded_off, ref_off);
1176 return nodded;
1177}
1178
1179/*----------------------------------------------------------------------------*/
1186/*----------------------------------------------------------------------------*/
1187static cpl_imagelist * sofi_spc_jitter_distor(
1188 cpl_imagelist * ilist,
1189 const char * arc)
1190{
1191 cpl_polynomial * arc_poly;
1192 cpl_polynomial * id_poly;
1193 cpl_table * tab;
1194 cpl_size power[2];
1195 cpl_vector * profile;
1196 cpl_imagelist * warped_list;
1197 cpl_image * warped;
1198 int i;
1199
1200 /* Test entries */
1201 if (ilist == NULL) return NULL;
1202 if (arc == NULL) return NULL;
1203
1204 /* Get the arc polynomial */
1205 arc_poly = cpl_polynomial_new(2);
1206 cpl_msg_info(cpl_func, "Get the arc distortion from the file");
1207 if ((tab = cpl_table_load(arc, 1, 0)) == NULL) {
1208 cpl_msg_error(cpl_func, "cannot load the arc table");
1209 cpl_polynomial_delete(arc_poly);
1210 return NULL;
1211 }
1212 for (i=0; i<cpl_table_get_nrow(tab); i++) {
1213 power[0] = cpl_table_get_int(tab, "Degree_of_x", i, NULL);
1214 power[1] = cpl_table_get_int(tab, "Degree_of_y", i, NULL);
1215 cpl_polynomial_set_coeff(arc_poly, power,
1216 cpl_table_get_double(tab, "poly2d_coef", i, NULL));
1217 }
1218 cpl_table_delete(tab);
1219
1220 /* Create the ID polynomial for the other direction */
1221 id_poly = cpl_polynomial_new(2);
1222 power[0] = 1;
1223 power[1] = 0;
1224 cpl_polynomial_set_coeff(id_poly, power, 1.0);
1225
1226 /* Create the kernel */
1227 profile = cpl_vector_new(CPL_KERNEL_DEF_SAMPLES);
1228 cpl_vector_fill_kernel_profile(profile, CPL_KERNEL_DEFAULT,
1229 CPL_KERNEL_DEF_WIDTH);
1230
1231 /* Correct the images */
1232 warped_list = cpl_imagelist_new();
1233 for (i=0; i<cpl_imagelist_get_size(ilist); i++) {
1234 warped = cpl_image_duplicate(cpl_imagelist_get(ilist, i));
1235 if (cpl_image_warp_polynomial(warped, cpl_imagelist_get(ilist, i),
1236 id_poly, arc_poly, profile, CPL_KERNEL_DEF_WIDTH, profile,
1237 CPL_KERNEL_DEF_WIDTH) != CPL_ERROR_NONE) {
1238 cpl_msg_error(cpl_func, "cannot correct the distortion");
1239 cpl_image_delete(warped);
1240 cpl_polynomial_delete(arc_poly);
1241 cpl_polynomial_delete(id_poly);
1242 cpl_vector_delete(profile);
1243 return NULL;
1244 }
1245 cpl_imagelist_set(warped_list, warped, i);
1246 }
1247 cpl_vector_delete(profile);
1248 cpl_polynomial_delete(arc_poly);
1249 cpl_polynomial_delete(id_poly);
1250 return warped_list;
1251}
1252
1253/*----------------------------------------------------------------------------*/
1260/*----------------------------------------------------------------------------*/
1261static double sofi_spc_jitter_refine_offset(
1262 cpl_image * ima1,
1263 cpl_image * ima2)
1264{
1265 double pos1, pos2;
1266
1267 /* Test entries */
1268 if (ima1 == NULL) return 10000.0;
1269 if (ima2 == NULL) return 10000.0;
1270
1271 /* Detect the spectra */
1272 if (irplib_spectrum_find_brightest(ima1, 0.0, NO_SHADOW, 0.0, 1,
1273 &pos1) == -1){
1274 return 10000.0;
1275 }
1276 if (irplib_spectrum_find_brightest(ima2, 0.0, NO_SHADOW, 0.0, 1,
1277 &pos2) == -1){
1278 return 10000.0;
1279 }
1280 return pos1-pos2;
1281}
1282
1283/*----------------------------------------------------------------------------*/
1289/*----------------------------------------------------------------------------*/
1290static cpl_table * sofi_spc_jitter_extract(cpl_image * combined)
1291{
1292 cpl_errorstate allstate = cpl_errorstate_get();
1293 int ri_dist, le_dist, ri_width, le_width, spec_pos;
1294 int nx, ny;
1295 double pos;
1296 int right_side, left_side;
1297 int sky_pos[4];
1298 cpl_vector * sky;
1299 cpl_vector * spec;
1300 cpl_vector * wl;
1301 double * pspec;
1302 double * psky;
1303 double * pwl;
1304 cpl_table * out;
1305 cpl_bivector * toplot;
1306 int throw;
1307 int res;
1308 int i;
1309 int nokfirst = -1;
1310 int noklast = 0;
1311
1312 /* Test entries */
1313 if (combined == NULL) return NULL;
1314
1315 /* Initialise */
1316 nx = cpl_image_get_size_x(combined);
1317 ny = cpl_image_get_size_y(combined);
1318 ri_dist = sofi_spc_jitter_config.extr_sky_ri_dist;
1319 le_dist = sofi_spc_jitter_config.extr_sky_le_dist;
1320 ri_width = sofi_spc_jitter_config.extr_sky_ri_width;
1321 le_width = sofi_spc_jitter_config.extr_sky_le_width;
1322 spec_pos = sofi_spc_jitter_config.extr_spec_pos;
1323 res = -1;
1324
1325 /* Detect the spectrum position if not passed */
1326 if (spec_pos < 0) {
1327 if (sofi_spc_jitter_config.throws == NULL) {
1328 cpl_msg_error(cpl_func,
1329 "Need a throw value to detect the spectra !!");
1330 return NULL;
1331 }
1332 for (i=0; i<cpl_vector_get_size(sofi_spc_jitter_config.throws); i++){
1333 throw = (int)cpl_vector_get(sofi_spc_jitter_config.throws, i);
1334 if ((res = irplib_spectrum_find_brightest(combined, throw,
1335 TWO_SHADOWS, 0.0, 1, &pos)) == 0) break;
1336 if ((res = irplib_spectrum_find_brightest(combined, throw,
1337 ONE_SHADOW, 0.0, 1, &pos)) == 0) break;
1338 }
1339 if (res != 0) {
1340 cpl_msg_error(cpl_func, "Cannot detect the spectrum");
1341 (void)cpl_error_set(cpl_func, CPL_ERROR_DATA_NOT_FOUND);
1342 return NULL;
1343 }
1344 spec_pos = (int)pos;
1345 cpl_msg_info(cpl_func, "Spectrum detected at x = %d", spec_pos);
1346 }
1347
1348 /* Set the parameters for the extraction */
1349
1350 /* Spectrum position */
1351 left_side = spec_pos - (int)(sofi_spc_jitter_config.extr_spec_width/2);
1352 right_side = left_side + sofi_spc_jitter_config.extr_spec_width;
1353 if ((left_side < 1) || (right_side > ny)) {
1354 cpl_msg_error(cpl_func, "Spectrum zone falls outside the image");
1355 (void)cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
1356 "%d < 1 || %d > %d", (int)left_side,
1357 (int)right_side, (int)ny);
1358 return NULL;
1359 }
1360 /* Residual Sky position */
1361 if (ri_dist < 0) ri_dist = 2*sofi_spc_jitter_config.extr_spec_width;
1362 if (le_dist < 0) le_dist = 2*sofi_spc_jitter_config.extr_spec_width;
1363 sky_pos[1] = spec_pos - le_dist;
1364 sky_pos[0] = sky_pos[1] - le_width;
1365 sky_pos[2] = spec_pos + ri_dist;
1366 sky_pos[3] = sky_pos[2] + ri_width;
1367
1368 /* Get the sky */
1369 sky = cpl_vector_new(ny);
1370 psky = cpl_vector_get_data(sky);
1371 if (((sky_pos[0] < 1) || (le_width == 0)) &&
1372 ((sky_pos[3] <= nx) && (ri_width > 0))) {
1373 for (i=0; i<ny; i++) {
1374 psky[i] = cpl_image_get_median_window(combined, sky_pos[2],
1375 ny-i, sky_pos[3], ny-i);
1376 }
1377 } else if (((sky_pos[3] > nx) || (ri_width == 0))
1378 && ((sky_pos[0] > 0) && (le_width > 0))) {
1379 for (i=0; i<ny; i++) {
1380 psky[i] = cpl_image_get_median_window(combined, sky_pos[0],
1381 ny-i, sky_pos[1], ny-i);
1382 }
1383 } else if ((ri_width != 0) && (le_width != 0)
1384 && (sky_pos[0] > 0) && (sky_pos[3] <= nx)) {
1385 for (i=0; i<ny; i++) {
1386 psky[i] = cpl_image_get_median_window(combined, sky_pos[2],
1387 ny-i, sky_pos[3], ny-i);
1388 psky[i] += cpl_image_get_median_window(combined, sky_pos[0],
1389 ny-i, sky_pos[1], ny-i);
1390 psky[i] /= 2.0;
1391 }
1392 } else {
1393 for (i=0; i<ny; i++) psky[i] = 0.0;
1394 }
1395
1396 /* Estimate the spectrum */
1397 spec = cpl_vector_new(ny);
1398 pspec = cpl_vector_get_data(spec);
1399 for (i=0; i<ny; i++) {
1400 cpl_errorstate prestate = cpl_errorstate_get();
1401 pspec[i] = cpl_image_get_flux_window(combined, left_side, ny-i,
1402 right_side, ny-i);
1403 if (cpl_errorstate_is_equal(prestate)) {
1404 if (nokfirst < 0) nokfirst = i;
1405 noklast = i;
1406 pspec[i] -= psky[i] * sofi_spc_jitter_config.extr_spec_width;
1407 } else {
1408 pspec[i] = 0.0;
1409 (void)cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
1410 "%d <= %d. i=%d < %d=ny", (int)left_side,
1411 (int)right_side, (int)i, (int)ny);
1412 }
1413 }
1414 if (nokfirst <= noklast) {
1415 cpl_msg_warning(cpl_func, "Only %d of %d column(s) processed, from %d "
1416 "to %d", (int)(1+noklast-nokfirst), (int)ny,
1417 (int)(ny-noklast), (int)(ny-nokfirst));
1418 cpl_errorstate_dump(allstate, CPL_FALSE,
1419 cpl_errorstate_dump_one_warning);
1420 cpl_errorstate_set(allstate);
1421 } else {
1422 (void)cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
1423 "All %d columns failed", (int)ny);
1424 }
1425
1426 /* Get the wavelength */
1427 wl = cpl_vector_new(ny);
1428 pwl = cpl_vector_get_data(wl);
1429 for (i=0; i<ny; i++) {
1430 pwl[i] = sofi_spc_jitter_config.wavecal_a0 +
1431 sofi_spc_jitter_config.wavecal_a1 * (i+1) +
1432 sofi_spc_jitter_config.wavecal_a2 * (i+1) * (i+1) +
1433 sofi_spc_jitter_config.wavecal_a3 * (i+1) * (i+1) * (i+1) +
1434 sofi_spc_jitter_config.wavecal_a4 * (i+1) * (i+1) * (i+1) + (i+1);
1435 }
1436
1437 /* Plot the spectrum if requested */
1438 if (sofi_spc_jitter_config.display) {
1439 toplot = cpl_bivector_wrap_vectors(wl, spec);
1440 cpl_plot_bivector(NULL, "t 'Spectrum' w lines", NULL, toplot);
1441 cpl_bivector_unwrap_vectors(toplot);
1442 toplot = cpl_bivector_wrap_vectors(wl, sky);
1443 cpl_plot_bivector(NULL, "t 'Sky' w lines", NULL, toplot);
1444 cpl_bivector_unwrap_vectors(toplot);
1445 }
1446
1447 /* Create and fill the output table */
1448 out = cpl_table_new(nx);
1449 cpl_table_new_column(out, "Y_coordinate", CPL_TYPE_DOUBLE);
1450 cpl_table_new_column(out, "Extracted_spectrum_value", CPL_TYPE_DOUBLE);
1451 cpl_table_new_column(out, "Sky_spectrum", CPL_TYPE_DOUBLE);
1452 for (i=0; i<nx; i++) {
1453 cpl_table_set_double(out, "Y_coordinate", i, pwl[i]);
1454 cpl_table_set_double(out, "Extracted_spectrum_value", i, pspec[i]);
1455 cpl_table_set_double(out, "Sky_spectrum", i, psky[i]);
1456 }
1457 cpl_vector_delete(wl);
1458 cpl_vector_delete(spec);
1459 cpl_vector_delete(sky);
1460 return out;
1461}
1462
1463/*----------------------------------------------------------------------------*/
1472/*----------------------------------------------------------------------------*/
1473static cpl_error_code sofi_spc_jitter_save(cpl_frameset * set,
1474 const cpl_image * ima,
1475 const cpl_table * tab,
1476 const cpl_parameterlist * parlist)
1477{
1478 cpl_propertylist * plist;
1479 /* Get the QC params in qclist */
1480 cpl_propertylist * qclist = cpl_propertylist_new();
1481 cpl_propertylist * paflist;
1482 /* Get the reference frame */
1483 const cpl_frame * ref_frame
1484 = irplib_frameset_get_first_from_group(set, CPL_FRAME_GROUP_RAW);
1485 char qc_str[128];
1486 int i;
1487
1488
1489 /* Add QC parameters */
1490 if (sofi_spc_jitter_config.filter[0] != (char)0)
1491 cpl_propertylist_append_string(qclist, "ESO QC FILTER OBS",
1492 sofi_spc_jitter_config.filter);
1493 cpl_propertylist_append_double(qclist, "ESO QC DISPCO1",
1494 sofi_spc_jitter_config.wavecal_a0);
1495 cpl_propertylist_append_double(qclist, "ESO QC DISPCO2",
1496 sofi_spc_jitter_config.wavecal_a1);
1497 cpl_propertylist_append_double(qclist, "ESO QC DISPCO3",
1498 sofi_spc_jitter_config.wavecal_a2);
1499 cpl_propertylist_append_double(qclist, "ESO QC DISPCO4",
1500 sofi_spc_jitter_config.wavecal_a3);
1501 cpl_propertylist_append_double(qclist, "ESO QC DISPCO5",
1502 sofi_spc_jitter_config.wavecal_a4);
1503 cpl_propertylist_append_double(qclist, "ESO QC WLEN",
1504 sofi_spc_jitter_config.wavecal_a0 +
1505 sofi_spc_jitter_config.wavecal_a1 * 512 +
1506 sofi_spc_jitter_config.wavecal_a2 * 512 * 512 +
1507 sofi_spc_jitter_config.wavecal_a3 * 512 * 512 * 512 +
1508 sofi_spc_jitter_config.wavecal_a4 * 512 * 512 * 512 * 512);
1509 cpl_propertylist_append_double(qclist, "ESO QC DISP XCORR",
1510 sofi_spc_jitter_config.wavecal_cc);
1511 if (sofi_spc_jitter_config.wavecal_out == 0) {
1512 cpl_propertylist_append_string(qclist, "ESO QC WLMETHOD",
1513 "physical model");
1514 } else if (sofi_spc_jitter_config.wavecal_out == 1) {
1515 cpl_propertylist_append_string(qclist, "ESO QC WLMETHOD",
1516 "sky lines");
1517 } else if (sofi_spc_jitter_config.wavecal_out == 2) {
1518 cpl_propertylist_append_string(qclist, "ESO QC WLMETHOD",
1519 "arc file");
1520 }
1521 for (i=0; i<cpl_vector_get_size(sofi_spc_jitter_config.intensities);i++) {
1522 sprintf(qc_str, "ESO QC SPEC INTENS%d", i+1);
1523 cpl_propertylist_append_double(qclist, qc_str,
1524 cpl_vector_get(sofi_spc_jitter_config.intensities, i));
1525 }
1526
1527 /* Change WCS keywords to the computed wavelength solution */
1528 cpl_propertylist_update_double(qclist, "CRVAL1",
1529 sofi_spc_jitter_config.wavecal_a0);
1530 cpl_propertylist_update_double(qclist, "CRVAL2", 1.0);
1531 cpl_propertylist_update_double(qclist, "CRPIX1", 1.0);
1532 cpl_propertylist_update_double(qclist, "CRPIX2", 1.0);
1533 cpl_propertylist_update_double(qclist, "CDELT1",
1534 sofi_spc_jitter_config.wavecal_a1);
1535 cpl_propertylist_update_double(qclist, "CDELT2", 1.0);
1536 cpl_propertylist_update_string(qclist, "CTYPE1", "LINEAR");
1537 cpl_propertylist_update_string(qclist, "CTYPE2", "LINEAR");
1538 if (cpl_propertylist_has(qclist, "CD1_1")) {
1539 cpl_propertylist_update_double(qclist, "CD1_1",
1540 sofi_spc_jitter_config.wavecal_a1);
1541 } else {
1542 cpl_propertylist_insert_after_double(qclist, "CTYPE2", "CD1_1",
1543 sofi_spc_jitter_config.wavecal_a1);
1544 }
1545 if (cpl_propertylist_has(qclist, "CD2_2")) {
1546 cpl_propertylist_update_double(qclist, "CD2_2", 1.0);
1547 } else {
1548 cpl_propertylist_insert_after_double(qclist, "CD1_1", "CD2_2", 1.0);
1549 }
1550
1551 /* Write the image */
1552 irplib_dfs_save_image(set, parlist, set, ima, CPL_BPP_IEEE_FLOAT,
1553 RECIPE_STRING, SOFI_SPC_JITTER_COMB, qclist,
1554 NULL, PACKAGE "/" PACKAGE_VERSION,
1555 RECIPE_STRING "_combined" CPL_DFS_FITS);
1556
1557 /* Write the table */
1558 if (tab != NULL) {
1559 irplib_dfs_save_table(set, parlist, set, tab, NULL, RECIPE_STRING,
1560 SOFI_SPC_JITTER_EXTR, qclist, NULL, PACKAGE "/" PACKAGE_VERSION,
1561 RECIPE_STRING "_extracted" CPL_DFS_FITS);
1562 }
1563
1564
1565 /* Get FITS header from reference file */
1566 plist = cpl_propertylist_load(cpl_frame_get_filename(ref_frame), 0);
1567 if (plist == NULL) {
1568 cpl_propertylist_delete(qclist);
1569 return cpl_error_set_message(cpl_func, cpl_error_get_code(), "Could "
1570 "not read header from reference frame");
1571 }
1572
1573 /* Get the keywords for the paf file */
1574 paflist = cpl_propertylist_new();
1575 cpl_propertylist_copy_property_regexp(paflist, plist,
1576 "^(ARCFILE|MJD-OBS|INSTRUME|ESO TPL ID|ESO TPL NEXP|ESO DPR CATG|"
1577 "ESO DPR TECH|ESO DPR TYPE|DATE-OBS|ESO INS GRAT NAME|"
1578 "ESO INS GRAT WLEN|ESO INS OPTI1 ID|ESO OBS ID|ESO OBS TARG NAME)$", 0);
1579 cpl_propertylist_delete(plist);
1580
1581 /* Copy the QC in paflist */
1582 cpl_propertylist_copy_property_regexp(paflist, qclist, ".", 0);
1583 cpl_propertylist_delete(qclist);
1584
1585 /* Save the PAF file */
1586 cpl_dfs_save_paf("SOFI", RECIPE_STRING, paflist,
1587 RECIPE_STRING CPL_DFS_PAF);
1588 cpl_propertylist_delete(paflist);
1589
1590 return cpl_error_get_code();
1591}
1592
1593/*----------------------------------------------------------------------------*/
1601/*----------------------------------------------------------------------------*/
1602static int off_comp(double off1, double off2, double thresh)
1603{
1604 if (((off1>thresh) && (off2<thresh)) || ((off1<thresh) && (off2>thresh)))
1605 return 1;
1606 else return 0;
1607}
1608
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_cumoffsetx
double sofi_pfits_get_cumoffsetx(const cpl_propertylist *plist)
find out the cumulative offset in X
Definition sofi_pfits.c:110
sofi_pfits_get_filter
const char * sofi_pfits_get_filter(const cpl_propertylist *plist)
find out which wave band is active in short wavelength
Definition sofi_pfits.c:243
sofi_extract_frameset
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
sofi_extract_filename
const char * sofi_extract_filename(const cpl_frameset *in, const char *tag)
Extract the filename ffor the first frame of the given tag.
Definition sofi_utils.c:342
sofi_correct_crosstalk_list
int sofi_correct_crosstalk_list(cpl_imagelist *ilist)
Remove the Cross-talk effect from an image list.
Definition sofi_utils.c:91
sofi_get_disprel_estimate
cpl_polynomial * sofi_get_disprel_estimate(const char *filename, int poly_deg)
Estimate the instrument dispersion relation.
Definition sofi_wavelength.c:313
sofi_get_slitwidth
double sofi_get_slitwidth(const char *filename)
Find out the slit width.
Definition sofi_wavelength.c:263
sofi_wavelength_engine
cpl_polynomial * sofi_wavelength_engine(const cpl_image *in, const char *table_name, const char *oh, const char *xe, const char *ne, const cpl_polynomial *phdisprel, double slit_width, int degree, double wl_err, int nsamples, int use_ppm, int plot, double *xcorr)
Compute a dispersion relation.
Definition sofi_wavelength.c:150
Generated on Thu Jun 18 2026 12:20:44 for IIINSTRUMENT Pipeline Reference Manual by doxygen 1.13.2