IIINSTRUMENT Pipeline Reference Manual 1.5.16
sofi_img_jitter.c
1/* $Id: sofi_img_jitter.c,v 1.30 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.30 $
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 <string.h>
38#include <cpl.h>
39
40#include "irplib_utils.h"
41#include "irplib_calib.h"
42
43#include "sofi_utils.h"
44#include "sofi_pfits.h"
45#include "sofi_dfs.h"
46
47/*-----------------------------------------------------------------------------
48 Define
49 -----------------------------------------------------------------------------*/
50
51#define NEGLIG_OFF_DIFF 0.1
52#define SQR(x) ((x)*(x))
53
54/*-----------------------------------------------------------------------------
55 Functions prototypes
56 -----------------------------------------------------------------------------*/
57
58static int sofi_img_jitter_create(cpl_plugin *);
59static int sofi_img_jitter_exec(cpl_plugin *);
60static int sofi_img_jitter_destroy(cpl_plugin *);
61static int sofi_img_jitter(cpl_parameterlist *, cpl_frameset *);
62
63static cpl_image ** sofi_img_jitter_reduce(cpl_frameset *, cpl_frameset *,
64 const char *, const char *, const char *, cpl_vector **);
65static cpl_imagelist ** sofi_img_jitter_load(cpl_frameset *, cpl_frameset *);
66static cpl_vector * sofi_img_jitter_sky(cpl_imagelist **, cpl_imagelist *);
67static cpl_vector * sofi_img_jitter_sky_running(cpl_imagelist **);
68static cpl_image ** sofi_img_jitter_saa(cpl_imagelist *,
69 cpl_frameset *);
70static int sofi_img_jitter_sub_row_median(cpl_image *);
71static cpl_table * sofi_img_jitter_qc(cpl_image *);
72static double sofi_img_jitter_get_mode(cpl_vector *);
73static int sofi_img_jitter_save(cpl_image *, cpl_table *, cpl_vector *,
74 cpl_parameterlist *, cpl_frameset *);
75
76/*-----------------------------------------------------------------------------
77 Static variables
78 -----------------------------------------------------------------------------*/
79
80static struct {
81 /* Inputs */
82 const char * offsets;
83 const char * objects;
84 int crosstalk;
85 int sky_minnb;
86 int sky_halfw;
87 int sky_rejmin;
88 int sky_rejmax;
89 int sx;
90 int sy;
91 int mx;
92 int my;
93 cpl_geom_combine comb_meth;
94 int rej_low;
95 int rej_high;
96 int row_med;
97 /* Outputs */
98 double pixscale;
99 double dit;
100 int nb_obj_frames;
101 int nb_sky_frames;
102 int nb_rej_frames;
103 double iq;
104 int nbobjs;
105 double fwhm_pix;
106 double fwhm_arcsec;
107 double fwhm_mode;
108} sofi_img_jitter_config;
109
110static char sofi_img_jitter_description[] =
111"sofi_img_jitter -- SOFI imaging jitter recipe.\n"
112"The files listed in the Set Of Frames (sof-file) must be tagged:\n"
113"raw-file.fits "SOFI_IMG_JITTER_OBJ_RAW" or\n"
114"raw-file.fits "SOFI_IMG_JITTER_SKY_RAW" or\n"
115"flat-file.fits "SOFI_CALIB_FLAT" or\n"
116"dark-file.fits "SOFI_CALIB_DARK"\n";
117
118/*-----------------------------------------------------------------------------
119 Functions code
120 -----------------------------------------------------------------------------*/
121
122/*----------------------------------------------------------------------------*/
130/*----------------------------------------------------------------------------*/
131int cpl_plugin_get_info(cpl_pluginlist * list)
132{
133 cpl_recipe * recipe = cpl_calloc(1, sizeof(*recipe));
134 cpl_plugin * plugin = &recipe->interface;
135
136 cpl_plugin_init(plugin,
137 CPL_PLUGIN_API,
138 SOFI_BINARY_VERSION,
139 CPL_PLUGIN_TYPE_RECIPE,
140 "sofi_img_jitter",
141 "Jitter recipe",
142 sofi_img_jitter_description,
143 "Yves Jung",
144 "yjung@eso.org",
145 sofi_get_license(),
146 sofi_img_jitter_create,
147 sofi_img_jitter_exec,
148 sofi_img_jitter_destroy);
149
150 cpl_pluginlist_append(list, plugin);
151
152 return 0;
153}
154
155/*----------------------------------------------------------------------------*/
164/*----------------------------------------------------------------------------*/
165static int sofi_img_jitter_create(cpl_plugin * plugin)
166{
167 cpl_recipe * recipe;
168 cpl_parameter * p;
169
170 /* Get the recipe out of the plugin */
171 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
172 recipe = (cpl_recipe *)plugin;
173 else return CPL_ERROR_UNSPECIFIED;
174
175 /* Create the parameters list in the cpl_recipe object */
176 recipe->parameters = cpl_parameterlist_new();
177
178 /* Fill the parameters list */
179 /* --off */
180 p = cpl_parameter_new_value("sofi.sofi_img_jitter.offsets", CPL_TYPE_STRING,
181 "Offsets file. If the i'th frame is believed "
182 "to have the offset (CumoffsetX_i, "
183 "CumoffsetY_i), then the i'th line in "
184 "the offset file should be:\n "
185 "X_i Y_i\n "
186 "where\n "
187 "X_i = CumoffsetX_0 - CumoffsetX_i\n "
188 "and\n "
189 "Y_i = CumoffsetY_0 - CumoffsetY_i.\n "
190 "(This implies that the first line in the "
191 "offset file should consist of two zeroes).",
192 "sofi.sofi_img_jitter", NULL);
193 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "off");
194 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
195 cpl_parameterlist_append(recipe->parameters, p);
196
197 /* --objs */
198 p = cpl_parameter_new_value("sofi.sofi_img_jitter.objects",
199 CPL_TYPE_STRING, "objects file", "sofi.sofi_img_jitter", NULL);
200 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "objs");
201 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
202 cpl_parameterlist_append(recipe->parameters, p);
203
204 /* --crosstalk */
205 p = cpl_parameter_new_value("sofi.sofi_img_jitter.crosstalk", CPL_TYPE_BOOL,
206 "flag to remove the crosstalk effect", "sofi.sofi_img_jitter",
207 TRUE);
208 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "crosstalk");
209 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
210 cpl_parameterlist_append(recipe->parameters, p);
211
212 /* --sky_par */
213 p = cpl_parameter_new_value("sofi.sofi_img_jitter.sky_par",
214 CPL_TYPE_STRING,
215 "Rejection parameters for sky filtering",
216 "sofi.sofi_img_jitter",
217 "10,7,3,3");
218 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sky_par");
219 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
220 cpl_parameterlist_append(recipe->parameters, p);
221
222 /* --xcorr */
223 p = cpl_parameter_new_value("sofi.sofi_img_jitter.xcorr",
224 CPL_TYPE_STRING,
225 "Cross correlation search and measure sizes",
226 "sofi.sofi_img_jitter",
227 "40,40,65,65");
228 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "xcorr");
229 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
230 cpl_parameterlist_append(recipe->parameters, p);
231
232 /* --comb_meth */
233 p = cpl_parameter_new_value("sofi.sofi_img_jitter.comb_meth",
234 CPL_TYPE_STRING, "union / inter / first", "sofi.sofi_img_jitter",
235 "union");
236 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "comb_meth");
237 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
238 cpl_parameterlist_append(recipe->parameters, p);
239
240 /* --rej */
241 p = cpl_parameter_new_value("sofi.sofi_img_jitter.rej",
242 CPL_TYPE_STRING,
243 "Low and high number of rejected values",
244 "sofi.sofi_img_jitter",
245 "2,2");
246 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "rej");
247 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
248 cpl_parameterlist_append(recipe->parameters, p);
249
250 /* --row_med */
251 p = cpl_parameter_new_value("sofi.sofi_img_jitter.row_med", CPL_TYPE_BOOL,
252 "flag to subtract the median of each row", "sofi.sofi_img_jitter",
253 TRUE);
254 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "row_med");
255 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
256 cpl_parameterlist_append(recipe->parameters, p);
257
258 /* Return */
259 return 0;
260}
261
262/*----------------------------------------------------------------------------*/
268/*----------------------------------------------------------------------------*/
269static int sofi_img_jitter_exec(cpl_plugin * plugin)
270{
271 cpl_recipe * recipe;
272
273 /* Get the recipe out of the plugin */
274 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
275 recipe = (cpl_recipe *)plugin;
276 else return CPL_ERROR_UNSPECIFIED;
277
278 return sofi_img_jitter(recipe->parameters, recipe->frames);
279}
280
281/*----------------------------------------------------------------------------*/
287/*----------------------------------------------------------------------------*/
288static int sofi_img_jitter_destroy(cpl_plugin * plugin)
289{
290 cpl_recipe * recipe;
291
292 /* Get the recipe out of the plugin */
293 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
294 recipe = (cpl_recipe *)plugin;
295 else return CPL_ERROR_UNSPECIFIED;
296
297 cpl_parameterlist_delete(recipe->parameters);
298 return 0;
299}
300
301/*----------------------------------------------------------------------------*/
308/*----------------------------------------------------------------------------*/
309static int sofi_img_jitter(
310 cpl_parameterlist * parlist,
311 cpl_frameset * framelist)
312{
313 const char * sval;
314 cpl_parameter * par;
315 const char * flat;
316 const char * dark;
317 cpl_frameset * objframes;
318 cpl_frameset * skyframes;
319 cpl_image ** combined;
320 cpl_table * objs_stats;
321 cpl_vector * sky_bg;
322
323 /* Initialise */
324 par = NULL;
325 sofi_img_jitter_config.pixscale = -1.0;
326 sofi_img_jitter_config.dit = -1.0;
327 sofi_img_jitter_config.iq = -1.0;
328 sofi_img_jitter_config.nbobjs = -1;
329 sofi_img_jitter_config.fwhm_pix = -1.0;
330 sofi_img_jitter_config.fwhm_arcsec = -1.0;
331 sofi_img_jitter_config.fwhm_mode = -1.0;
332 sofi_img_jitter_config.offsets = NULL;
333 sofi_img_jitter_config.objects = NULL;
334 sofi_img_jitter_config.nb_obj_frames = 0;
335 sofi_img_jitter_config.nb_rej_frames = 0;
336 sofi_img_jitter_config.nb_sky_frames = 0;
337
338 /* Retrieve input parameters */
339 /* Offsets */
340 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.offsets");
341 sofi_img_jitter_config.offsets = cpl_parameter_get_string(par);
342 /* Objects */
343 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.objects");
344 sofi_img_jitter_config.objects = cpl_parameter_get_string(par);
345 /* Cross-talk */
346 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.crosstalk");
347 sofi_img_jitter_config.crosstalk = cpl_parameter_get_bool(par);
348 /* Rejection parameters for sky filtering */
349 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.sky_par");
350 sval = cpl_parameter_get_string(par);
351 if (sscanf(sval, "%d,%d,%d,%d",
352 &sofi_img_jitter_config.sky_minnb,
353 &sofi_img_jitter_config.sky_halfw,
354 &sofi_img_jitter_config.sky_rejmin,
355 &sofi_img_jitter_config.sky_rejmax)!=4) {
356 return CPL_ERROR_UNSPECIFIED;
357 }
358 /* Cross correlation windows parameters */
359 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.xcorr");
360 sval = cpl_parameter_get_string(par);
361 if (sscanf(sval, "%d,%d,%d,%d",
362 &sofi_img_jitter_config.sx,
363 &sofi_img_jitter_config.sy,
364 &sofi_img_jitter_config.mx,
365 &sofi_img_jitter_config.my)!=4) {
366 return CPL_ERROR_UNSPECIFIED;
367 }
368 /* --comb_meth */
369 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.comb_meth");
370 sval = cpl_parameter_get_string(par);
371 if (!strcmp(sval, "union"))
372 sofi_img_jitter_config.comb_meth = CPL_GEOM_UNION;
373 else if (!strcmp(sval, "inter"))
374 sofi_img_jitter_config.comb_meth = CPL_GEOM_INTERSECT;
375 else if (!strcmp(sval, "first"))
376 sofi_img_jitter_config.comb_meth = CPL_GEOM_FIRST;
377 else {
378 cpl_msg_error(cpl_func, "Invalid combine method specified");
379 return CPL_ERROR_UNSPECIFIED;
380 }
381
382 /* Number of rejected values in stacking */
383 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.rej");
384 sval = cpl_parameter_get_string(par);
385 if (sscanf(sval, "%d,%d",
386 &sofi_img_jitter_config.rej_low,
387 &sofi_img_jitter_config.rej_high)!=2) {
388 return CPL_ERROR_UNSPECIFIED;
389 }
390 /* Row median */
391 par = cpl_parameterlist_find(parlist, "sofi.sofi_img_jitter.row_med");
392 sofi_img_jitter_config.row_med = cpl_parameter_get_bool(par);
393
394 /* Identify the RAW and CALIB frames in the input frameset */
395 if (sofi_dfs_set_groups(framelist)) {
396 cpl_msg_error(cpl_func, "Cannot identify RAW and CALIB frames");
397 return CPL_ERROR_UNSPECIFIED;
398 }
399
400 /* Retrieve calibration data */
401 flat = sofi_extract_filename(framelist, SOFI_CALIB_FLAT);
402 dark = sofi_extract_filename(framelist, SOFI_CALIB_DARK);
403
404 /* Retrieve raw frames */
405 if ((objframes = sofi_extract_frameset(framelist,
406 SOFI_IMG_JITTER_OBJ_RAW)) == NULL) {
407 cpl_msg_error(cpl_func, "Cannot find objs frames in the input list");
408 return CPL_ERROR_UNSPECIFIED;
409 }
410 skyframes = sofi_extract_frameset(framelist, SOFI_IMG_JITTER_SKY_RAW);
411
412 /* Apply the reduction */
413 cpl_msg_info(cpl_func, "Apply the data recombination");
414 cpl_msg_indent_more();
415 if ((combined = sofi_img_jitter_reduce(objframes, skyframes, flat,
416 dark, NULL, &sky_bg)) == NULL) {
417 cpl_msg_error(cpl_func, "Cannot recombine the data");
418 cpl_frameset_delete(objframes);
419 if (skyframes) cpl_frameset_delete(skyframes);
420 cpl_msg_indent_less();
421 return CPL_ERROR_UNSPECIFIED;
422 }
423 cpl_frameset_delete(objframes);
424 if (skyframes) cpl_frameset_delete(skyframes);
425 cpl_msg_indent_less();
426
427 /* Compute QC parameters from the combined image */
428 cpl_msg_info(cpl_func, "Compute QC parameters from the combined image");
429 cpl_msg_indent_more();
430 if ((objs_stats = sofi_img_jitter_qc(combined[0])) == NULL) {
431 cpl_msg_warning(cpl_func, "Cannot compute all parameters");
432 }
433 cpl_msg_indent_less();
434
435 /* Save the products */
436 cpl_msg_info(cpl_func, "Save the products");
437 cpl_msg_indent_more();
438 if (sofi_img_jitter_save(combined[0], objs_stats, sky_bg, parlist,
439 framelist) == -1){
440 cpl_msg_error(cpl_func, "Cannot save the products");
441 cpl_image_delete(combined[0]);
442 cpl_image_delete(combined[1]);
443 cpl_free(combined);
444 if (objs_stats) cpl_table_delete(objs_stats);
445 if (sky_bg) cpl_vector_delete(sky_bg);
446 cpl_msg_indent_less();
447 return CPL_ERROR_UNSPECIFIED;
448 }
449 cpl_msg_indent_less();
450
451 /* Return */
452 cpl_image_delete(combined[0]);
453 cpl_image_delete(combined[1]);
454 cpl_free(combined);
455 if (objs_stats) cpl_table_delete(objs_stats);
456 if (sky_bg) cpl_vector_delete(sky_bg);
457 return 0;
458}
459
460/*----------------------------------------------------------------------------*/
471/*----------------------------------------------------------------------------*/
472static cpl_image ** sofi_img_jitter_reduce(
473 cpl_frameset * obj,
474 cpl_frameset * sky,
475 const char * flat,
476 const char * dark,
477 const char * bpm,
478 cpl_vector ** skybg)
479{
480 cpl_imagelist ** in;
481 cpl_image ** combined;
482
483 /* Initialise */
484 *skybg = NULL;
485
486 /* Load the input data */
487 cpl_msg_info(cpl_func, "Load the input data");
488 cpl_msg_indent_more();
489 if ((in = sofi_img_jitter_load(obj, sky)) == NULL) {
490 cpl_msg_error(cpl_func, "Cannot load input data");
491 cpl_msg_indent_less();
492 return NULL;
493 }
494 cpl_msg_indent_less();
495
496 /* Apply the cross-talk correction */
497 if (sofi_img_jitter_config.crosstalk) {
498 cpl_msg_info(cpl_func, "Apply the cross-talk correction");
499 cpl_msg_indent_more();
500 if (sofi_correct_crosstalk_list(in[0]) == -1) {
501 cpl_msg_error(cpl_func, "Cannot correct for Cross-talk");
502 cpl_imagelist_delete(in[0]);
503 if (in[1]) cpl_imagelist_delete(in[1]);
504 cpl_free(in);
505 cpl_msg_indent_less();
506 return NULL;
507 }
508 if (in[1]) {
509 if (sofi_correct_crosstalk_list(in[1]) == -1) {
510 cpl_msg_error(cpl_func, "Cannot correct for Cross-talk");
511 cpl_imagelist_delete(in[0]);
512 if (in[1]) cpl_imagelist_delete(in[1]);
513 cpl_free(in);
514 cpl_msg_indent_less();
515 return NULL;
516 }
517 }
518 cpl_msg_indent_less();
519 }
520
521 /* Apply the calibrations */
522 if (flat || dark || bpm) {
523 cpl_msg_info(cpl_func, "Apply the calibrations");
524 cpl_msg_indent_more();
525 if (irplib_flat_dark_bpm_calib(in[0], flat, dark, bpm) == -1) {
526 /* Calibrate the objects */
527 cpl_msg_error(cpl_func, "Cannot calibrate the objects");
528 cpl_imagelist_delete(in[0]);
529 if (in[1]) cpl_imagelist_delete(in[1]);
530 cpl_free(in);
531 cpl_msg_indent_less();
532 return NULL;
533 }
534 if (in[1]) {
535 if (irplib_flat_dark_bpm_calib(in[1], flat, dark, bpm) == -1) {
536 /* Calibrate the sky */
537 cpl_msg_error(cpl_func, "Cannot calibrate the sky");
538 cpl_imagelist_delete(in[0]);
539 if (in[1]) cpl_imagelist_delete(in[1]);
540 cpl_free(in);
541 cpl_msg_indent_less();
542 return NULL;
543 }
544 }
545 cpl_msg_indent_less();
546 }
547
548 /* Apply the sky correction */
549 cpl_msg_info(cpl_func, "Sky estimation and correction");
550 cpl_msg_indent_more();
551 if ((*skybg = sofi_img_jitter_sky(&(in[0]), in[1])) == NULL) {
552 cpl_msg_error(cpl_func, "Cannot estimate the sky");
553 cpl_imagelist_delete(in[0]);
554 if (in[1]) cpl_imagelist_delete(in[1]);
555 cpl_free(in);
556 cpl_msg_indent_less();
557 return NULL;
558 }
559 cpl_msg_indent_less();
560 sofi_img_jitter_config.nb_obj_frames = cpl_imagelist_get_size(in[0]);
561 if (in[1] != NULL)
562 sofi_img_jitter_config.nb_sky_frames = cpl_imagelist_get_size(in[1]);
563 if (in[1]) cpl_imagelist_delete(in[1]);
564 in[1] = NULL;
565
566 /* Apply the shift and add */
567 cpl_msg_info(cpl_func, "Shift and add");
568 cpl_msg_indent_more();
569 combined = sofi_img_jitter_saa(in[0], obj);
570 if (combined == NULL) {
571 cpl_msg_error(cpl_func, "Cannot apply the shift and add");
572 cpl_imagelist_delete(in[0]);
573 cpl_free(in);
574 if (*skybg != NULL) cpl_vector_delete(*skybg);
575 *skybg = NULL;
576 cpl_msg_indent_less();
577 return NULL;
578 }
579 cpl_imagelist_delete(in[0]);
580 cpl_free(in);
581 cpl_msg_indent_less();
582
583 /* Post processing on the combined image */
584 if (sofi_img_jitter_config.row_med) {
585 cpl_msg_info(cpl_func, "Subtract the median from each row");
586 sofi_img_jitter_sub_row_median(combined[0]);
587 }
588
589 return combined;
590}
591
592/*----------------------------------------------------------------------------*/
599/*----------------------------------------------------------------------------*/
600static cpl_imagelist ** sofi_img_jitter_load(
601 cpl_frameset * obj,
602 cpl_frameset * sky)
603{
604 cpl_imagelist ** isets;
605 cpl_frame * frame;
606 cpl_propertylist * plist;
607
608 /* Test entries */
609 if (obj == NULL) return NULL;
610
611 /* Check if ok */
612 if (cpl_error_get_code()) return NULL;
613
614 /* Get the frame type pixscale and DIT value */
615 frame = cpl_frameset_get_position(obj, 0);
616 plist=cpl_propertylist_load(cpl_frame_get_filename(frame), 0);
617 sofi_img_jitter_config.pixscale = sofi_pfits_get_pixscale(plist);
618 sofi_img_jitter_config.dit = sofi_pfits_get_dit(plist);
619 cpl_propertylist_delete(plist);
620 if (cpl_error_get_code()) {
621 cpl_msg_error(cpl_func, "Missing keyword in FITS header");
622 return NULL;
623 }
624
625 /* Allocate the image sets */
626 isets = cpl_malloc(2 * sizeof(cpl_imagelist *));
627 isets[0] = cpl_imagelist_load_frameset(obj, CPL_TYPE_FLOAT, 1, 0);
628 if (sky != NULL)
629 isets[1] = cpl_imagelist_load_frameset(sky, CPL_TYPE_FLOAT, 1, 0);
630 else isets[1] = NULL;
631
632 /* Check if at least the objects are there */
633 if (isets[0] == NULL) {
634 cpl_msg_error(cpl_func, "The objects frames could not be loaded");
635 if (isets[1] != NULL) cpl_imagelist_delete(isets[1]);
636 cpl_free(isets);
637 return NULL;
638 }
639 /* Return */
640 return isets;
641}
642
643/*----------------------------------------------------------------------------*/
650/*----------------------------------------------------------------------------*/
651static cpl_vector * sofi_img_jitter_sky(
652 cpl_imagelist ** objs,
653 cpl_imagelist * skys)
654{
655 cpl_image * sky;
656 cpl_vector * bg;
657 int sky_method;
658 int nframes, nskys;
659 double median;
660 cpl_image * cur_ima;
661 int i;
662
663 /* Initialise */
664 nframes = cpl_imagelist_get_size(*objs);
665 bg = NULL;
666
667 /* Decide the sky method */
668 if (sofi_img_jitter_config.sky_minnb > nframes) sky_method = 1;
669 else sky_method = 2;
670
671 /* Compute the sky frame */
672 if (skys != NULL) {
673 cpl_msg_info(cpl_func, "Median of sky images");
674 /* Get the median of the sky images */
675 nskys = cpl_imagelist_get_size(skys);
676 bg = cpl_vector_new(nskys);
677 for (i=0 ; i<nskys ; i++) {
678 median = cpl_image_get_median(cpl_imagelist_get(skys, i));
679 cpl_vector_set(bg, i, median);
680 }
681 /* Use sky images */
682 if ((sky = cpl_imagelist_collapse_median_create(skys)) == NULL) {
683 cpl_msg_error(cpl_func, "Cannot compute the median of sky images");
684 cpl_vector_delete(bg);
685 return NULL;
686 }
687 /* Correct the objects images */
688 if (cpl_imagelist_subtract_image(*objs, sky) != CPL_ERROR_NONE) {
689 cpl_msg_error(cpl_func, "Cannot correct the obj images from the sky");
690 cpl_image_delete(sky);
691 cpl_vector_delete(bg);
692 return NULL;
693 }
694 cpl_image_delete(sky);
695 /* Normalise the object planes */
696 for (i=0 ; i<nframes ; i++) {
697 cur_ima = cpl_imagelist_get(*objs, i);
698 median = cpl_image_get_median(cur_ima);
699 cpl_image_subtract_scalar(cur_ima, median);
700 }
701 } else if (sky_method == 1) {
702 cpl_msg_info(cpl_func, "Median of object images");
703 /* Use objs images */
704 if ((sky = cpl_imagelist_collapse_median_create(*objs)) == NULL) {
705 cpl_msg_error(cpl_func, "Cannot compute the median of obj images");
706 return NULL;
707 }
708 /* Correct the objects images */
709 if (cpl_imagelist_subtract_image(*objs, sky) != CPL_ERROR_NONE) {
710 cpl_msg_error(cpl_func, "Cannot correct the obj images from the sky");
711 cpl_image_delete(sky);
712 return NULL;
713 }
714 /* Create a 1 value vector */
715 bg = cpl_vector_new(1);
716 cpl_vector_set(bg, 0, cpl_image_get_median(sky));
717 cpl_image_delete(sky);
718 /* Normalise the object planes */
719 for (i=0 ; i<nframes ; i++) {
720 cur_ima = cpl_imagelist_get(*objs, i);
721 median = cpl_image_get_median(cur_ima);
722 cpl_image_subtract_scalar(cur_ima, median);
723 }
724 } else if (sky_method == 2) {
725 cpl_msg_info(cpl_func, "Running filter on object images");
726 /* Use objects images */
727 if ((bg = sofi_img_jitter_sky_running(objs)) == NULL) {
728 cpl_msg_error(cpl_func, "Cannot apply the running filter for the sky");
729 return NULL;
730 }
731 }
732 /* Free and return */
733 return bg;
734}
735
736/*----------------------------------------------------------------------------*/
755/*----------------------------------------------------------------------------*/
756static cpl_vector * sofi_img_jitter_sky_running(cpl_imagelist ** in)
757{
758 int rejmin, rejmax, halfw;
759 cpl_imagelist * filtres;
760 int ni, nx, ny, pos;
761 cpl_vector * medians;
762 cpl_image * cur_ima;
763 float * pcur_ima;
764 cpl_image * tmp_ima;
765 float * ptmp_ima;
766 cpl_vector * localwin;
767 int fr_p, to_p, n_curp;
768 cpl_vector * bg;
769 double bg_val, out;
770 float one_med;
771 int i, j, k;
772
773 /* Test entries */
774 if (in==NULL || *in == NULL) return NULL;
775
776 /* Initialise */
777 rejmin = sofi_img_jitter_config.sky_rejmin;
778 rejmax = sofi_img_jitter_config.sky_rejmax;
779 halfw = sofi_img_jitter_config.sky_halfw;
780 ni = cpl_imagelist_get_size(*in);
781 cur_ima = cpl_imagelist_get(*in, 0);
782 nx = cpl_image_get_size_x(cur_ima);
783 ny = cpl_image_get_size_y(cur_ima);
784
785 /* Tests on validity of rejection parameters */
786 if (((rejmin+rejmax)>=halfw) || (halfw<1) || (rejmin<0) || (rejmax<0)) {
787 cpl_msg_error(cpl_func, "cannot run filter with rej parms %d (%d-%d)",
788 halfw, rejmin, rejmax);
789 return NULL;
790 }
791 /* Pre-compute median value in each plane */
792 medians = cpl_vector_new(ni);
793 for (i=0 ; i<ni ; i++) {
794 cur_ima = cpl_imagelist_get(*in, i);
795 cpl_vector_set(medians, i, cpl_image_get_median(cur_ima));
796 }
797 /* Allocate output cube */
798 filtres = cpl_imagelist_new();
799
800 /* Allocate output bg */
801 bg = cpl_vector_new(ni);
802
803 /* Main loop over input planes */
804 for (k=0 ; k<ni ; k++) {
805 /* Allocate output plane */
806 tmp_ima = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
807 ptmp_ima = cpl_image_get_data_float(tmp_ima);
808 /* Compute border indices */
809 fr_p = k - halfw;
810 to_p = k + halfw;
811 if (fr_p<0) fr_p=0;
812 if (to_p>(ni-1)) to_p=ni-1;
813
814 /* Number of valid planes to consider after edge effects */
815 n_curp = to_p - fr_p;
816
817 /* Allocate local window */
818 localwin = cpl_vector_new(n_curp);
819
820 bg_val = 0.0;
821 /* Loop over all pixels */
822 for (pos=0 ; pos<nx*ny ; pos++) {
823 /* Fill up local window */
824 j=0;
825 for (i=fr_p ; i<=to_p ; i++) {
826 if (i!=k) {
827 cur_ima = cpl_imagelist_get(*in, i);
828 pcur_ima = cpl_image_get_data_float(cur_ima);
829 cpl_vector_set(localwin, j,
830 (double)pcur_ima[pos]-cpl_vector_get(medians, i));
831 j++;
832 }
833 }
834 /* Sort window */
835 cpl_vector_sort(localwin, 1);
836 /* Reject min and max, accumulate other pixels */
837 out = 0.0;
838 for (i=rejmin ; i<(n_curp-rejmax) ; i++) {
839 out += cpl_vector_get(localwin, i);
840 }
841 /* Take the mean */
842 out /= (double)(n_curp - rejmin - rejmax);
843 /* Assign value */
844 cur_ima = cpl_imagelist_get(*in, k);
845 pcur_ima = cpl_image_get_data_float(cur_ima);
846 ptmp_ima[pos] = pcur_ima[pos] -
847 (float)(out + cpl_vector_get(medians, k));
848
849 bg_val += (out+cpl_vector_get(medians, k));
850 }
851 cpl_vector_delete(localwin);
852 cpl_vector_set(bg, k, bg_val/(nx*ny));
853 cpl_imagelist_set(filtres, tmp_ima, k);
854 }
855 cpl_imagelist_delete(*in);
856 cpl_vector_delete(medians);
857
858 /* Subtract median from each frame */
859 for (i=0 ; i<ni ; i++) {
860 cur_ima = cpl_imagelist_get(filtres, i);
861 one_med = cpl_image_get_median(cur_ima);
862 cpl_image_subtract_scalar(cur_ima, one_med);
863 }
864 *in = filtres;
865 return bg;
866}
867
868/*----------------------------------------------------------------------------*/
875/*----------------------------------------------------------------------------*/
876static cpl_image ** sofi_img_jitter_saa(
877 cpl_imagelist * in,
878 cpl_frameset * objframes)
879{
880 cpl_frame * frame;
881 cpl_propertylist * plist;
882 cpl_bivector * offsets_est;
883 double * offsets_est_x;
884 double * offsets_est_y;
885 cpl_bivector * objs;
886 double * objs_x;
887 double * objs_y;
888 cpl_apertures * aperts;
889 cpl_image ** combined;
890 cpl_vector * thresh_vect;
891 cpl_image * diff;
892 int nfiles;
893 int i;
894
895 /* Get the number of images */
896 nfiles = cpl_imagelist_get_size(in);
897 if (cpl_frameset_get_size(objframes) != nfiles) {
898 cpl_msg_error(cpl_func, "Invalid input objects sizes");
899 return NULL;
900 }
901
902 /* Get the offsets estimation of each input file pair */
903 cpl_msg_info(cpl_func, "Get the offsets estimation");
904 offsets_est = NULL;
905 if (sofi_img_jitter_config.offsets &&
906 sofi_img_jitter_config.offsets[0] != (char)0) {
907 /* A file has been provided on the command line */
908 offsets_est = cpl_bivector_read(sofi_img_jitter_config.offsets);
909 if ((offsets_est==NULL)||(cpl_bivector_get_size(offsets_est)!=nfiles)) {
910 cpl_msg_error(cpl_func, "Cannot get offsets from %s",
911 sofi_img_jitter_config.offsets);
912 return NULL;
913 }
914 } else {
915 /* Get the offsets from the header */
916 offsets_est = cpl_bivector_new(nfiles);
917 offsets_est_x = cpl_bivector_get_x_data(offsets_est);
918 offsets_est_y = cpl_bivector_get_y_data(offsets_est);
919 for (i=0 ; i<nfiles ; i++) {
920 if (cpl_error_get_code()) {
921 cpl_bivector_delete(offsets_est);
922 cpl_ensure(0, cpl_error_get_code(), NULL);
923 }
924 /* X and Y offsets */
925 frame = cpl_frameset_get_position(objframes, i);
926 plist=cpl_propertylist_load(cpl_frame_get_filename(frame),0);
927 offsets_est_x[i] = -1.0 * sofi_pfits_get_cumoffsetx(plist);
928 offsets_est_y[i] = -1.0 * sofi_pfits_get_cumoffsety(plist);
929 cpl_propertylist_delete(plist);
930 if (cpl_error_get_code()) {
931 cpl_msg_warning(cpl_func, "Cannot get offsets from header");
932 cpl_bivector_delete(offsets_est);
933 cpl_ensure(0, cpl_error_get_code(), NULL);
934 }
935 }
936 /* Subtract the first offset to all offsets */
937 for (i=1 ; i<nfiles ; i++) {
938 offsets_est_x[i] -= offsets_est_x[0];
939 offsets_est_y[i] -= offsets_est_y[0];
940 }
941 offsets_est_x[0] = offsets_est_y[0] = 0.00;
942 }
943
944 /* Read the provided objects file if provided */
945 objs = NULL;
946 if (sofi_img_jitter_config.objects &&
947 sofi_img_jitter_config.objects[0] != (char)0) {
948 cpl_msg_info(cpl_func, "Get the user provided correlation objects");
949 /* A file has been provided on the command line */
950 objs = cpl_bivector_read(sofi_img_jitter_config.objects);
951 if (objs==NULL) {
952 cpl_msg_error(cpl_func, "Cannot get objects from %s",
953 sofi_img_jitter_config.objects);
954 if (offsets_est) cpl_bivector_delete(offsets_est);
955 return NULL;
956 }
957 }
958
959 /* Get a correlation point from the difference of the first images */
960 if (objs == NULL) {
961 cpl_msg_info(cpl_func, "Get a cross-correlation point");
962 thresh_vect = cpl_vector_new(4);
963 cpl_vector_set(thresh_vect, 0, 5.0);
964 cpl_vector_set(thresh_vect, 1, 2.0);
965 cpl_vector_set(thresh_vect, 2, 1.0);
966 cpl_vector_set(thresh_vect, 3, 0.5);
967 diff = cpl_image_subtract_create(cpl_imagelist_get(in, 0),
968 cpl_imagelist_get(in, 1));
969 if ((aperts = cpl_apertures_extract_window(diff, thresh_vect,
970 200, 200, 800, 800, NULL)) == NULL) {
971 cpl_msg_error(cpl_func, "Cannot find any cross-correlation point");
972 if (offsets_est) cpl_bivector_delete(offsets_est);
973 cpl_vector_delete(thresh_vect);
974 cpl_image_delete(diff);
975 return NULL;
976 }
977 cpl_image_delete(diff);
978 cpl_vector_delete(thresh_vect);
979 cpl_apertures_sort_by_npix(aperts);
980 objs = cpl_bivector_new(1);
981 objs_x = cpl_bivector_get_x_data(objs);
982 objs_y = cpl_bivector_get_y_data(objs);
983 objs_x[0] = cpl_apertures_get_pos_x(aperts, 1);
984 objs_y[0] = cpl_apertures_get_pos_y(aperts, 1);
985 cpl_apertures_delete(aperts);
986 if (objs == NULL) {
987 cpl_msg_error(cpl_func, "Cannot find any cross-correlation point");
988 if (offsets_est) cpl_bivector_delete(offsets_est);
989 return NULL;
990 }
991 cpl_msg_info(cpl_func, "Correlation point: %g %g\n", objs_x[0], objs_y[0]);
992 }
993
994 /* Recombine the images */
995 cpl_msg_info(cpl_func, "Recombine the images set");
996 cpl_msg_indent_more();
997 if ((combined = cpl_geom_img_offset_combine(in, offsets_est, 1, objs,
998 NULL, NULL,
999 sofi_img_jitter_config.sx,
1000 sofi_img_jitter_config.sy,
1001 sofi_img_jitter_config.mx,
1002 sofi_img_jitter_config.my,
1003 sofi_img_jitter_config.rej_low,
1004 sofi_img_jitter_config.rej_high,
1005 sofi_img_jitter_config.comb_meth)) == NULL) {
1006 cpl_msg_error(cpl_func, "Cannot recombine the images");
1007 if (offsets_est) cpl_bivector_delete(offsets_est);
1008 cpl_bivector_delete(objs);
1009 cpl_msg_indent_less();
1010 return NULL;
1011 }
1012 /* Update QC params */
1013 i = (int)(cpl_image_get_max(combined[1]));
1014 sofi_img_jitter_config.nb_rej_frames =
1015 sofi_img_jitter_config.nb_obj_frames - i;
1016 sofi_img_jitter_config.nb_obj_frames = i;
1017 cpl_msg_indent_less();
1018
1019 /* Free and return */
1020 if (offsets_est) cpl_bivector_delete(offsets_est);
1021 cpl_bivector_delete(objs);
1022 return combined;
1023}
1024
1025/*----------------------------------------------------------------------------*/
1031/*----------------------------------------------------------------------------*/
1032static int sofi_img_jitter_sub_row_median(cpl_image * in)
1033{
1034 float * pin;
1035 int nx, ny;
1036 double median;
1037 int i, j;
1038
1039 /* Test entries */
1040 if (in == NULL) return CPL_ERROR_UNSPECIFIED;
1041
1042 /* Initialise */
1043 nx = cpl_image_get_size_x(in);
1044 ny = cpl_image_get_size_y(in);
1045 pin = cpl_image_get_data_float(in);
1046
1047 for (i=0 ; i<ny ; i++) {
1048 median = cpl_image_get_median_window(in, 1, i+1, nx, i+1);
1049 for (j=0 ; j<nx ; j++) {
1050 pin[j+i*nx] -= (float)median;
1051 }
1052 }
1053 return 0;
1054}
1055
1056/*----------------------------------------------------------------------------*/
1062/*----------------------------------------------------------------------------*/
1063static cpl_table * sofi_img_jitter_qc(cpl_image * combined)
1064{
1065 cpl_apertures * aperts;
1066 cpl_bivector * fwhms;
1067 double * fwhms_x;
1068 double * fwhms_y;
1069 cpl_vector * fwhms_good;
1070 double * fwhms_good_data;
1071 int nb_val, nb_good;
1072 double f_min, f_max, fr, fx, fy;
1073 cpl_vector * thresh_vec;
1074 cpl_table * out_tab;
1075 int nb_objs;
1076 int i, j;
1077
1078 /* Initialise */
1079 double seeing_min_arcsec = 0.1;
1080 double seeing_max_arcsec = 5.0;
1081 double seeing_fwhm_var = 0.2;
1082
1083 /* Check entries */
1084 if (combined == NULL) return NULL;
1085
1086 /* Create the vector for the detection thresholds */
1087 thresh_vec = cpl_vector_new(4);
1088 cpl_vector_set(thresh_vec, 0, 5.0);
1089 cpl_vector_set(thresh_vec, 1, 2.0);
1090 cpl_vector_set(thresh_vec, 2, 1.0);
1091 cpl_vector_set(thresh_vec, 3, 0.5);
1092
1093 /* Detect apertures */
1094 if ((aperts = cpl_apertures_extract(combined, thresh_vec, NULL)) == NULL) {
1095 cpl_msg_error(cpl_func, "Cannot detect any aperture");
1096 cpl_vector_delete(thresh_vec);
1097 return NULL;
1098 }
1099 cpl_vector_delete(thresh_vec);
1100
1101 /* Number of detected objects */
1102 nb_objs = cpl_apertures_get_size(aperts);
1103 sofi_img_jitter_config.nbobjs = nb_objs;
1104
1105 /* Compute the FHWM of the detected apertures */
1106 CPL_DIAG_PRAGMA_PUSH_IGN(-Wdeprecated-declarations);
1107 fwhms = cpl_apertures_get_fwhm(combined, aperts);
1108 CPL_DIAG_PRAGMA_POP;
1109
1110 if (fwhms == NULL) {
1111 cpl_msg_error(cpl_func, "Cannot compute the FWHMs");
1112 cpl_apertures_delete(aperts);
1113 return NULL;
1114 }
1115
1116 /* Access the data */
1117 nb_val = cpl_bivector_get_size(fwhms);
1118 fwhms_x = cpl_bivector_get_x_data(fwhms);
1119 fwhms_y = cpl_bivector_get_y_data(fwhms);
1120
1121 /* Store everything in the output table */
1122 out_tab = cpl_table_new(nb_objs);
1123 cpl_table_new_column(out_tab, "POS_X", CPL_TYPE_DOUBLE);
1124 cpl_table_new_column(out_tab, "POS_Y", CPL_TYPE_DOUBLE);
1125 cpl_table_new_column(out_tab, "FWHM_X", CPL_TYPE_DOUBLE);
1126 cpl_table_new_column(out_tab, "FWHM_Y", CPL_TYPE_DOUBLE);
1127 cpl_table_new_column(out_tab, "FLUX", CPL_TYPE_DOUBLE);
1128 for (i=0 ; i<nb_objs ; i++) {
1129 cpl_table_set_double(out_tab, "POS_X", i,
1130 cpl_apertures_get_centroid_x(aperts, i+1));
1131 cpl_table_set_double(out_tab, "POS_Y", i,
1132 cpl_apertures_get_centroid_y(aperts, i+1));
1133 cpl_table_set_double(out_tab, "FWHM_X", i, fwhms_x[i]);
1134 cpl_table_set_double(out_tab, "FWHM_Y", i, fwhms_y[i]);
1135 cpl_table_set_double(out_tab, "FLUX", i,
1136 cpl_apertures_get_flux(aperts, i+1));
1137 }
1138 cpl_apertures_delete(aperts);
1139
1140 /* Get the number of good values */
1141 nb_good = 0;
1142 for (i=0 ; i<nb_val ; i++) {
1143 if ((fwhms_x[i] > 0.0) && (fwhms_y[i] > 0.0)) nb_good++;
1144 }
1145 if (nb_good == 0) {
1146 cpl_table_delete(out_tab);
1147 cpl_bivector_delete(fwhms);
1148 (void)cpl_error_set(cpl_func, CPL_ERROR_DATA_NOT_FOUND);
1149 return NULL;
1150 }
1151
1152 /* Get the good values */
1153 fwhms_good = cpl_vector_new(nb_good);
1154 fwhms_good_data = cpl_vector_get_data(fwhms_good);
1155 j=0;
1156 for (i=0 ; i<nb_val ; i++) {
1157 if ((fwhms_x[i] > 0.0) && (fwhms_y[i] > 0.0)) {
1158 fwhms_good_data[j] = (fwhms_x[i]+fwhms_y[i])/2.0;
1159 j++;
1160 }
1161 }
1162
1163 /* Compute the fwhm */
1164 if (nb_good < 3) {
1165 /* Too few values to compute the median */
1166 sofi_img_jitter_config.fwhm_pix = fwhms_good_data[0];
1167 } else {
1168 /* Compute the median */
1169 sofi_img_jitter_config.fwhm_pix = cpl_vector_get_median_const(fwhms_good);
1170 }
1171 sofi_img_jitter_config.fwhm_arcsec =
1172 sofi_img_jitter_config.fwhm_pix * sofi_img_jitter_config.pixscale;
1173
1174 /* Compute the mode of the FWHMs */
1175 if (nb_good > 5) {
1176 sofi_img_jitter_config.fwhm_mode =
1177 sofi_img_jitter_get_mode(fwhms_good);
1178 }
1179 cpl_vector_delete(fwhms_good);
1180
1181 /* IQ is the median of the (fwhm_x+fwhm_y/2) of the good stars */
1182 /* Compute f_min and f_max */
1183 f_min = seeing_min_arcsec / sofi_img_jitter_config.pixscale;
1184 f_max = seeing_max_arcsec / sofi_img_jitter_config.pixscale;
1185
1186 /* Get the number of good values */
1187 nb_good = 0;
1188 for (i=0 ; i<nb_val ; i++) {
1189 fx = fwhms_x[i];
1190 fy = fwhms_y[i];
1191 fr = 2.0 * fabs(fx-fy) / (fx+fy);
1192 if ((fx > f_min) && (fx < f_max) && (fy > f_min) && (fy < f_max) &&
1193 (fr < seeing_fwhm_var)) nb_good++;
1194 }
1195 if (nb_good == 0) {
1196 cpl_table_delete(out_tab);
1197 cpl_bivector_delete(fwhms);
1198 (void)cpl_error_set(cpl_func, CPL_ERROR_DATA_NOT_FOUND);
1199 return NULL;
1200 }
1201
1202 /* Get the good values */
1203 fwhms_good = cpl_vector_new(nb_good);
1204 fwhms_good_data = cpl_vector_get_data(fwhms_good);
1205 j=0;
1206 for (i=0 ; i<nb_val ; i++) {
1207 fx = fwhms_x[i];
1208 fy = fwhms_y[i];
1209 fr = 2.0 * fabs(fx-fy) / (fx+fy);
1210 if ((fx > f_min) && (fx < f_max) && (fy > f_min) && (fy < f_max) &&
1211 (fr < seeing_fwhm_var)) {
1212 fwhms_good_data[j] = (fx + fy)/2.0;
1213 j++;
1214 }
1215 }
1216 cpl_bivector_delete(fwhms);
1217
1218 /* Compute the fwhm */
1219 if (nb_good < 3) {
1220 /* Too few values to compute the median */
1221 sofi_img_jitter_config.iq = fwhms_good_data[0];
1222 } else {
1223 /* Compute the median */
1224 sofi_img_jitter_config.iq = cpl_vector_get_median_const(fwhms_good);
1225 }
1226 cpl_vector_delete(fwhms_good);
1227 sofi_img_jitter_config.iq *= sofi_img_jitter_config.pixscale;
1228
1229 return out_tab;
1230}
1231
1232/*----------------------------------------------------------------------------*/
1238/*----------------------------------------------------------------------------*/
1239static double sofi_img_jitter_get_mode(cpl_vector * vec)
1240{
1241 int nb;
1242 int nbins;
1243 double min, max;
1244 double bin_size;
1245 cpl_bivector * hist;
1246 cpl_vector * hist_x;
1247 cpl_vector * hist_y;
1248 double cur_val;
1249 int cur_bin;
1250 double max_val;
1251 int max_bin;
1252 double mode;
1253 int i;
1254
1255 /* Test entries */
1256 if (vec == NULL) return -1.0;
1257
1258 /* Initialise */
1259 nb = cpl_vector_get_size(vec);
1260
1261 /* Create the histogram */
1262 nbins = 10;
1263 min = cpl_vector_get_min(vec);
1264 max = cpl_vector_get_max(vec);
1265 bin_size = (max-min)/nbins;
1266 hist = cpl_bivector_new(nbins);
1267 hist_x = cpl_bivector_get_x(hist);
1268 hist_y = cpl_bivector_get_y(hist);
1269 cpl_vector_fill(hist_x, 0.0);
1270 cpl_vector_fill(hist_y, 0.0);
1271 for (i=0 ; i<nbins ; i++) {
1272 cpl_vector_set(hist_x, i, min + i * bin_size);
1273 }
1274 for (i=0 ; i<nb ; i++) {
1275 cur_val = cpl_vector_get(vec, i);
1276 cur_bin = (int)((cur_val - min) / bin_size);
1277 if (cur_bin >= nbins) cur_bin -= 1.0;
1278 cur_val = cpl_vector_get(hist_y, cur_bin);
1279 cur_val += 1.0;
1280 cpl_vector_set(hist_y, cur_bin, cur_val);
1281 }
1282
1283 /* Get the mode of the histogram */
1284 max_val = cpl_vector_get(hist_y, 0);
1285 max_bin = 0;
1286 for (i=0 ; i<nbins ; i++) {
1287 cur_val = cpl_vector_get(hist_y, i);
1288 if (cur_val > max_val) {
1289 max_val = cur_val;
1290 max_bin = i;
1291 }
1292 }
1293 mode = cpl_vector_get(hist_x, max_bin);
1294 cpl_bivector_delete(hist);
1295 return mode;
1296}
1297
1298/*----------------------------------------------------------------------------*/
1308/*----------------------------------------------------------------------------*/
1309static int sofi_img_jitter_save(
1310 cpl_image * combined,
1311 cpl_table * objs_stats,
1312 cpl_vector * sky_bg,
1313 cpl_parameterlist * parlist,
1314 cpl_frameset * set)
1315{
1316 const cpl_frame * ref_frame;
1317 cpl_propertylist * plist;
1318 cpl_propertylist * paflist;
1319 cpl_propertylist * qclist;
1320 double pscale, dit, bg_mean, bg_stdev, bg_instmag;
1321 cpl_table * sky_bg_tab;
1322 int i;
1323
1324 /* Get the QC params in qclist */
1325 qclist = cpl_propertylist_new();
1326
1327 if (sky_bg != NULL) {
1328 pscale = sofi_img_jitter_config.pixscale;
1329 dit = sofi_img_jitter_config.dit;
1330 bg_mean = cpl_vector_get_mean(sky_bg);
1331 if (cpl_vector_get_size(sky_bg) < 2) bg_stdev = 0;
1332 else bg_stdev = cpl_vector_get_stdev(sky_bg);
1333 bg_instmag = -2.5 * log(bg_mean/(pscale*pscale*dit));
1334 cpl_propertylist_append_double(qclist, "ESO QC BACKGD MEAN", bg_mean);
1335 cpl_propertylist_append_double(qclist, "ESO QC BACKGD STDEV", bg_stdev);
1336 cpl_propertylist_append_double(qclist, "ESO QC BACKGD INSTMAG",
1337 bg_instmag);
1338 }
1339 cpl_propertylist_append_int(qclist, "ESO QC NBOBJS",
1340 sofi_img_jitter_config.nbobjs);
1341 cpl_propertylist_append_double(qclist, "ESO QC IQ",
1342 sofi_img_jitter_config.iq);
1343 cpl_propertylist_append_double(qclist, "ESO QC FWHM PIX",
1344 sofi_img_jitter_config.fwhm_pix);
1345 cpl_propertylist_append_double(qclist, "ESO QC FWHM ARCSEC",
1346 sofi_img_jitter_config.fwhm_arcsec);
1347 cpl_propertylist_append_double(qclist, "ESO QC FWHM MODE",
1348 sofi_img_jitter_config.fwhm_mode);
1349 cpl_propertylist_append_int(qclist, "ESO QC NB_OBJ_F",
1350 sofi_img_jitter_config.nb_obj_frames);
1351 cpl_propertylist_append_int(qclist, "ESO QC NB_SKY_F",
1352 sofi_img_jitter_config.nb_sky_frames);
1353 cpl_propertylist_append_int(qclist, "ESO QC NB_REJ_F",
1354 sofi_img_jitter_config.nb_rej_frames);
1355
1356 /* Write the image */
1357 irplib_dfs_save_image(set,
1358 parlist,
1359 set,
1360 combined,
1361 CPL_BPP_IEEE_FLOAT,
1362 "sofi_img_jitter",
1363 SOFI_IMG_JITTER_COMB,
1364 qclist,
1365 NULL,
1366 PACKAGE "/" PACKAGE_VERSION,
1367 "sofi_img_jitter.fits");
1368
1369
1370 /* Write the FITS background table */
1371 if (sky_bg) {
1372 /* First create the table from the vector */
1373 sky_bg_tab = cpl_table_new(cpl_vector_get_size(sky_bg));
1374 cpl_table_new_column(sky_bg_tab, "SKY_BG", CPL_TYPE_DOUBLE);
1375 for (i=0 ; i<cpl_vector_get_size(sky_bg) ; i++) {
1376 cpl_table_set_double(sky_bg_tab, "SKY_BG", i,
1377 cpl_vector_get(sky_bg, i));
1378 }
1379 irplib_dfs_save_table(set,
1380 parlist,
1381 set,
1382 sky_bg_tab,
1383 NULL,
1384 "sofi_img_jitter",
1385 SOFI_IMG_JITTER_BG,
1386 qclist,
1387 NULL,
1388 PACKAGE "/" PACKAGE_VERSION,
1389 "sofi_img_jitter_bg.fits");
1390 cpl_table_delete(sky_bg_tab);
1391 }
1392
1393 /* Write the FITS table */
1394 if (objs_stats) {
1395 irplib_dfs_save_table(set,
1396 parlist,
1397 set,
1398 objs_stats,
1399 NULL,
1400 "sofi_img_jitter",
1401 SOFI_IMG_JITTER_STARS,
1402 qclist,
1403 NULL,
1404 PACKAGE "/" PACKAGE_VERSION,
1405 "sofi_img_jitter_stars.fits");
1406 }
1407
1408 /* Get the reference frame */
1409 ref_frame = irplib_frameset_get_first_from_group(set, CPL_FRAME_GROUP_RAW);
1410
1411 /* Get FITS header from reference file */
1412 if ((plist=cpl_propertylist_load(cpl_frame_get_filename(ref_frame),
1413 0)) == NULL) {
1414 cpl_msg_error(cpl_func, "getting header from reference frame");
1415 cpl_propertylist_delete(qclist);
1416 return CPL_ERROR_UNSPECIFIED;
1417 }
1418
1419 /* Get the keywords for the paf file */
1420 paflist = cpl_propertylist_new();
1421 cpl_propertylist_copy_property_regexp(paflist, plist,
1422 "^(ARCFILE|MJD-OBS|INSTRUME|ESO TPL ID|ESO TPL NEXP|ESO DPR CATG|"
1423 "ESO DPR TECH|ESO DPR TYPE|DATE-OBS|ESO OBS ID|ESO INS PIXSCALE)$", 0);
1424 cpl_propertylist_delete(plist);
1425
1426 /* Copy the QC in paflist */
1427 cpl_propertylist_copy_property_regexp(paflist, qclist, ".", 0);
1428 cpl_propertylist_delete(qclist);
1429
1430 /* Save the PAF file */
1431 cpl_dfs_save_paf("SOFI",
1432 "sofi_img_jitter",
1433 paflist,
1434 "sofi_img_jitter.paf");
1435 cpl_propertylist_delete(paflist);
1436 return 0;
1437}
1438
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_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_cumoffsety
double sofi_pfits_get_cumoffsety(const cpl_propertylist *plist)
find out the cumulative offset in Y
Definition sofi_pfits.c:122
sofi_pfits_get_pixscale
double sofi_pfits_get_pixscale(const cpl_propertylist *plist)
find out the pixel scale
Definition sofi_pfits.c:449
sofi_get_license
const char * sofi_get_license(void)
Get the pipeline copyright and license.
Definition sofi_utils.c:60
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
Generated on Thu Jun 18 2026 12:20:44 for IIINSTRUMENT Pipeline Reference Manual by doxygen 1.13.2