/*---------------------------------------------------------------------------- File name : jitter_engine.c Author : N. Devillard Created on : Sept 30, 1997 Description : Jitter mode data reduction engine for ISAAC/SOFI ---------------------------------------------------------------------------*/ /* $Id: jitter_engine.c,v 1.76 2002/03/11 10:02:36 ndevilla Exp $ $Author: ndevilla $ $Date: 2002/03/11 10:02:36 $ $Revision: 1.76 $ */ /*---------------------------------------------------------------------------- Includes ---------------------------------------------------------------------------*/ #include #include "jitter_engine.h" /*---------------------------------------------------------------------------- Defines ---------------------------------------------------------------------------*/ #define XCORR_SRC_ERROR -1 #define XCORR_SRC_NOACQ 0 #define XCORR_SRC_AUTO 1 #define XCORR_SRC_FILE 2 #define X_ABS(x) (((x)>0)?(x):(-(x))) #define IBOOL(i) ((i)? "yes" : "no") #define IFTYP(i) (((i)==FRAME_OBJ) ? "obj" : "sky") /*---------------------------------------------------------------------------- Function codes ---------------------------------------------------------------------------*/ static int jitter_get_object_source(cube_t *, jitter_config *) ; static int jitter_remove_sky_background(cube_t **, jitter_config*); static int jitter_offset_skysub(cube_t **, jitter_config *); static int jitter_load_objects_from_file(jitter_config *); static int jitter_autofindobjects(cube_t *, jitter_config *); static int jitter_output_object_list(jitter_config *); static int jitter_find_offsets(cube_t **, jitter_config *); static int jitter_shiftandadd(cube_t **, jitter_config *); static int jitter_postproc(cube_t *, jitter_config *) ; static void jitter_add_comments(jitter_config *, qfits_header * fh); static int jitter_qc(jitter_config *, cube_t *) ; static int jitter_save_cube(cube_t * jit_cube, jitter_config *); /*---------------------------------------------------------------------------*/ /** @name jitter_engine @memo Launch the jitter engine from a given configuration. @param ini_name Name of ini file to use to initialize the engine @return int 0 if Ok, anything else otherwise @doc Main jitter engine function. The whole algorithm can be parsed here. */ /*---------------------------------------------------------------------------*/ #define NPARTS 8 long jitter_engine(char * ini_name) { jitter_config * cfg ; cube_t * jit_cube ; time_t local_t ; long total_pixels_in ; char status_file[FILENAMESZ] ; char filename[FILENAMESZ] ; int p ; p=0 ; /* Parse config file */ e_comment(0, "parsing configuration file...") ; if ((cfg = parse_jitter_ini_file(ini_name)) == NULL) { e_error("in parsing %s: aborting jitter", ini_name) ; return -1 ; } e_comment(0, "---> STARTING JITTER ENGINE") ; time(&local_t) ; e_comment(0, "%s", ctime(&local_t)) ; e_comment(0, "pid is %ld", (long)getpid()); /* * LOAD INPUT CUBE */ p++ ; e_comment(0, "---> part %d of %d: loading data", p, NPARTS) ; /* Test if there are enough frames */ if (cfg->nframes < MIN_NFRAMES_JITTER) { e_error("Not enough frames: %d<%d", cfg->nframes, MIN_NFRAMES_JITTER) ; cfg->status_load = ALGO_FAILED ; jitter_cfg_destroy(cfg) ; return 1 ; } /* Check that input frames come from ISAAC */ if (isaac_check_instrument(cfg->frames_filelist) == -1) { e_error("Data is not from ISAAC: [%s]", cfg->frames_filelist) ; cfg->status_load = ALGO_FAILED ; jitter_cfg_destroy(cfg) ; return 1 ; } /* Load the data */ if ((jit_cube = cube_load(cfg->frames_filelist)) == NULL) { e_error("loading input list: aborting") ; cfg->status_load = ALGO_FAILED ; jitter_cfg_destroy(cfg) ; return 1 ; } cfg->status_load = ALGO_OK ; /* Update number of pixels to be processed */ total_pixels_in = (long)(jit_cube->lx * jit_cube->ly * jit_cube->np) ; /* * APPLY DARK SUBTRACTION/FLAT-FIELD DIVISION */ p++ ; e_comment(0, "---> part %d of %d: dark/flat-field/bad pixels", p, NPARTS) ; if ((cfg->dark_subtraction == 0) && (cfg->flatfield_division == 0) && (cfg->badpixels_replacement == 0)) { cfg->status_cosmetics = ALGO_SKIPPED ; } else { isaac_ff_dark_badpix_handling(&jit_cube, cfg->flatfield_filename, cfg->dark_filename, cfg->badpixels_filename) ; cfg->status_cosmetics = ALGO_OK ; } /* * APPLY SKY BACKGROUND SUBTRACTION */ p++ ; e_comment(0, "---> part %d of %d: sky estimation/subtraction", p, NPARTS) ; if (jitter_remove_sky_background(&jit_cube, cfg)!=0) { e_error("in sky background subtraction: aborting") ; if (jit_cube!=NULL) cube_del(jit_cube); jitter_cfg_destroy(cfg); return -1 ; } /* * Find out objects in first frame, for x-correlation */ p++ ; e_comment(0, "---> part %d of %d: object acquisition", p, NPARTS); if (jitter_get_object_source(jit_cube, cfg)!=0) { e_error("in object acquisition: aborting"); cube_del(jit_cube); jitter_cfg_destroy(cfg); return -1 ; } /* * Find out the offsets */ p++ ; e_comment(0, "---> part %d of %d: offset determination", p, NPARTS) ; if (jitter_find_offsets(&jit_cube, cfg)!=0) { e_error("in offset finding section: aborting") ; cube_del(jit_cube) ; jitter_cfg_destroy(cfg); return -1 ; } /* Apply shift and add */ p++ ; e_comment(0, "---> part %d of %d: shift and add", p, NPARTS); if (jitter_shiftandadd(&jit_cube, cfg)!=0) { e_error("in shift and add"); cube_del(jit_cube); jitter_cfg_destroy(cfg); return -1 ; } /* * Optional post-processing features will be inserted here */ p++ ; e_comment(0, "---> part %d of %d: optional post-processing", p, NPARTS) ; if (jitter_postproc(jit_cube, cfg)!=0) { e_error("in post-processing") ; } /* * Save results */ p++ ; e_comment(0, "---> part %d of %d: saving output data", p, NPARTS); jitter_save_cube(jit_cube, cfg) ; cube_del(jit_cube) ; e_comment(0, "---> STOPPING JITTER ENGINE") ; time(&local_t) ; e_comment(0, "%s", ctime(&local_t)); /* If requested: launch an image viewer on the result */ if ((cfg->postproc_startviewer == 1) && (cfg->postproc_active == 1)) { sprintf(filename, "%s.fits", cfg->output_basename); show_image(filename, cfg->postproc_viewer) ; } /* If requested, output a status file as basename_status.ascii */ if (cfg->postproc_statusreport) { sprintf(status_file, "%s_status.ascii", cfg->output_basename); jitter_cfg_dump(cfg, status_file); } jitter_cfg_destroy(cfg); return total_pixels_in ; } /*-------------------------------------------------------------------------*/ /** @name jitter_remove_sky_background @memo Estimate and subtract the sky in the jitter cube. @param jit_cube Jitter cube. @param cfg Jitter blackboard. @return int 0 if Ok, anything else otherwise. @doc This function estimates and subtracts the sky in the jitter cube. */ /*--------------------------------------------------------------------------*/ static int jitter_remove_sky_background( cube_t ** jit_cube, jitter_config * cfg) { cube_t * object ; cube_t * skysub ; image_t * single_sky ; int i ; char output_name[1024] ; int status ; qfits_header* fh ; /* If no sky estimation is requested, bail out */ if (cfg->sky_estimate==0) { e_comment(1, "skipping sky estimation/subtraction") ; cfg->status_skyengine = ALGO_SKIPPED ; cfg->status_skybg = ALGO_SKIPPED ; return 0 ; } /* Sky estimation is done through a running filter */ /* * Special case: jitter+offset. * The input cube contains both sky and object frames. Only * one sky subtraction method is supported in this case: * median average the sky frames and subtract this average * frame from all objects. * The sky background is not computed in this case. */ if (cfg->contains_sky == 1) { cfg->status_skybg = ALGO_SKIPPED ; e_comment(1, "applying jitter+offset sky subtraction"); status = jitter_offset_skysub(jit_cube, cfg); if (status==0) { cfg->status_skyengine = ALGO_OK ; } else { cfg->status_skyengine = ALGO_FAILED ; } return status ; } object = *jit_cube ; /* * Check out if there are enough frames to perform sky * estimation If not, the estimated sky will just be a median * average of the input frames. */ if (object->np < cfg->sky_minnumofframes) { e_warning("not enough frames to use 3d filter: using median"); single_sky = cube_avg_median(object) ; if (single_sky == NULL) { e_error("cannot median average: aborting sky estimation") ; cfg->status_skyengine = ALGO_FAILED ; return 1 ; } skysub = cube_new(object->lx, object->ly, object->np); for (i=0 ; inp ; i++) { skysub->plane[i] = image_sub(object->plane[i], single_sky); } image_del(single_sky) ; cube_del(object) ; *jit_cube = skysub ; if (cfg->sky_outputdiff) { e_comment(1, "saving sky-subtracted frames"); sprintf(output_name, "%s_dif.fits", cfg->output_basename); /* Read the FITS header of the input file */ fh = qfits_header_read(cfg->input_list[0]) ; isaac_header_for_image(fh) ; isaac_pro_fits(fh, output_name, "REDUCED", NULL, isaac_imag_sw_jitter_diff, "OK", "img_jitter", cfg->nframes, NULL, NULL) ; cube_save_fits_hdrdump(skysub, output_name, fh) ; qfits_header_destroy(fh) ; e_comment(1, "Difference produced: [%s]", output_name) ; } cfg->status_skyengine = ALGO_OK ; cfg->status_skybg = ALGO_SKIPPED ; return 0 ; } /* * Normal sky combination */ cfg->sky_background = malloc((*jit_cube)->np * sizeof(double)); if (cfg->sky_sepquad) { e_comment(1, "separating quadrants for sky subtraction"); status = cube_3dfilt_runminmax_by_quad(jit_cube, cfg->sky_rejhalfwidth, cfg->sky_rejmin, cfg->sky_rejmax, cfg->sky_background) ; } else { status = cube_3dfilt_runminmax(jit_cube, cfg->sky_rejhalfwidth, cfg->sky_rejmin, cfg->sky_rejmax, cfg->sky_background) ; } if (status != 0) { e_error("in sky estimation: aborting") ; cfg->status_skyengine = ALGO_FAILED ; return status ; } /* Optional difference output */ if (cfg->sky_outputdiff) { e_comment(1, "saving sky-subtracted frames"); sprintf(output_name, "%s_dif.fits", cfg->output_basename); /* Read the FITS header of the input file */ fh = qfits_header_read(cfg->input_list[0]) ; isaac_header_for_image(fh) ; isaac_pro_fits(fh, output_name, "REDUCED", NULL, isaac_imag_sw_jitter_diff, "OK", "img_jitter", cfg->nframes, NULL, NULL) ; cube_save_fits_hdrdump(*jit_cube, output_name, fh) ; qfits_header_destroy(fh) ; e_comment(1, "Difference produced: [%s]", output_name) ; } cfg->status_skyengine = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_get_object_source @memo Find where x-correlation objects are coming from. @param jit_cube Input cube @param cfg Jitter blackboard. @return 0 if ok. @doc Check out that object detection is really needed, i.e. shift and add has been requested, and the provided offsets need refining or they are automatically computed. Object detection is also not needed if the objects are user provided. */ /*--------------------------------------------------------------------------*/ static int jitter_get_object_source( cube_t * jit_cube, jitter_config * cfg) { int need_acq ; int status ; if (cfg->saa_apply==0) { e_comment(1, "no shift and add requested: skipping object acq"); cfg->status_objacq = ALGO_SKIPPED ; return 0 ; } if ((cfg->saa_offinput == OFFSIN_HEADER)|| (cfg->saa_offinput == OFFSIN_FILE)) { need_acq = cfg->saa_off_refine ; } else { need_acq = 1 ; } if (need_acq == 0) { e_comment(1, "no refining requested: no need to detect objects"); cfg->status_objacq = ALGO_SKIPPED ; return 0 ; } if (cfg->saa_objacq_source==OBJACQ_FILE) { e_comment(1, "loading objects from user source"); status = jitter_load_objects_from_file(cfg); if (status==0) { cfg->status_objacq = ALGO_OK ; jitter_output_object_list(cfg); } else { cfg->status_objacq = ALGO_FAILED ; } return status ; } e_comment(1, "automatic object acquisition for x-correlation"); status = jitter_autofindobjects(jit_cube, cfg); if (status==0) { cfg->status_objacq = ALGO_OK ; jitter_output_object_list(cfg); } else { cfg->status_objacq = ALGO_FAILED ; } return status ; } /*-------------------------------------------------------------------------*/ /** @name jitter_autofindobjects @memo Find automatically relevant objects for x-correlation. @param jit_cube Jitter cube. @param cfg Jitter blackboard. @return int 0 if Ok, -1 otherwise. @doc This function find objects suitable for x-correlation in the requested detection image. The detected objects are placed in the blackboard. */ /*--------------------------------------------------------------------------*/ static int jitter_autofindobjects(cube_t * jit_cube, jitter_config * cfg) { image_t * detect_image ; image_t * second_image ; double3 * peaks ; int i ; int first_obj, first_sky ; /* * Associate detect_image with the correct frame */ switch (cfg->saa_autodetectimage) { case DETIM_FIRST: detect_image = image_load(cfg->input_list[0]); break ; case DETIM_DIFF: if (cfg->contains_sky==0) { detect_image = image_load(cfg->input_list[0]); second_image = image_load(cfg->input_list[1]); image_sub_local(detect_image, second_image); image_del(second_image); } else { first_obj = -1 ; first_sky = -1 ; i=0 ; while ((first_obj<0) || (first_sky<0)) { if (cfg->frametype[i] == FRAME_OBJ) first_obj=i ; if (cfg->frametype[i] == FRAME_SKY) first_sky=i ; if (i>(cfg->nframes-1)) { e_error("internal: in autofindobjects"); return -1 ; } i++; } if (first_obj<0) { e_error("cannot find any object in input??"); return -1 ; } if (first_sky<0) { e_error("cannot find any sky in input??"); return -1 ; } detect_image = image_load(cfg->input_list[first_obj]); second_image = image_load(cfg->input_list[first_sky]); image_sub_local(detect_image, second_image); image_del(second_image); } break ; default: e_error("internal: unknown detection image label"); return -1 ; } if (detect_image==NULL) { e_error("cannot get image for detection: aborting"); return -1 ; } /* * Detect suitable objects for x-correlation. */ peaks = get_xcorrelation_points(detect_image, cfg->saa_off_hx + cfg->saa_off_sx, cfg->saa_off_hy + cfg->saa_off_sy, cfg->saa_autothreshold, cfg->saa_autominpts, cfg->saa_automaxpts); image_del(detect_image); if (peaks==NULL) { e_error("cannot find enough valid points for xcorrelation"); return -1 ; } /* * Update relevant fields in config table */ cfg->xcorr_nobj = peaks->n ; cfg->xcorrobj_x = malloc(peaks->n * sizeof(int)) ; cfg->xcorrobj_y = malloc(peaks->n * sizeof(int)) ; for (i=0 ; in ; i++) { cfg->xcorrobj_x[i] = (int)peaks->x[i] ; cfg->xcorrobj_y[i] = (int)peaks->y[i] ; } double3_del(peaks); return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_load_objects_from_file @memo Load a list of objects from an ASCII file. @param cfg Jitter blackboard. @return int 0 if Ok, -1 otherwise. @doc This function loads a list of objects suitable for x-correlation from an ASCII file. The returned objects are stored in the blackboard. */ /*--------------------------------------------------------------------------*/ static int jitter_load_objects_from_file(jitter_config * cfg) { char * obj_list_name ; double3 * obj_list ; int i ; /* * Get the name of the object list from the cfg table */ obj_list_name = cfg->saa_objacq_filename ; e_comment(1, "user object file: [%s]", obj_list_name); if (!file_exists(obj_list_name)) { e_error("cannot find User Object list [%s]", obj_list_name); return -1 ; } /* Read the file using double3_read() */ obj_list = double3_read(obj_list_name); if (obj_list==NULL) { e_error("cannot read object list from [%s]", obj_list_name); return -1 ; } /* Save data into blackboard */ cfg->xcorr_nobj = obj_list->n ; cfg->xcorrobj_x = malloc(obj_list->n*sizeof(int)); cfg->xcorrobj_y = malloc(obj_list->n*sizeof(int)); for (i=0 ; in ; i++) { cfg->xcorrobj_x[i] = (int)(obj_list->x[i]+0.5) ; cfg->xcorrobj_y[i] = (int)(obj_list->y[i]+0.5) ; } /* Deallocate double3 list */ double3_del(obj_list); return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_output_object_list @memo Output the list of x-correlation objects to a file. @param cfg Jitter blackboard. @return int 0 if Ok, anything else otherwise. @doc Print out a list of x-correlation objects to an ASCII file. Does nothing if this option is not activated in the blackboard. */ /*--------------------------------------------------------------------------*/ static int jitter_output_object_list(jitter_config * cfg) { char obj_file_name[FILENAMESZ] ; FILE * obj_file ; int i ; char * s ; if (cfg->saa_outputobjs == 0) return 0 ; sprintf(obj_file_name, "%s_obj.paf", cfg->output_basename); if (cfg->xcorr_nobj<1) { e_error("strange number of objects [%d]: aborting output", cfg->xcorr_nobj); return -1 ; } if ((obj_file = paf_print_header(obj_file_name, "ISAAC/jitter", "List of objects")) == NULL) { e_error("cannot create objects PAF file") ; return -1 ; } /* Add ARCFILE */ if ((s = isaac_get_arcfile(cfg->input_list[0])) != NULL) fprintf(obj_file, "ARCFILE \"%s\" ;#\n", s) ; /* Add PRO.CATG */ fprintf(obj_file, "PRO.CATG \"%s\" ;# Product category\n", isaac_get_pro_catg_value(isaac_imag_sw_jitter_objs)) ; /* Add the date if available */ s = isaac_get_date_obs(cfg->input_list[0]) ; if (s) { fprintf(obj_file, "DATE-OBS \"%s\" ;# Date\n", s); } /* Add the filter ID */ s = isaac_get_filter(cfg->input_list[0]); if (s) { fprintf(obj_file, "INS.FILTER.ID \"%s\" ;# Filter name\n", s); } /* Add the Obs ID */ s = isaac_get_obs_id(cfg->input_list[0]); if (s) { fprintf(obj_file, "OBS.ID \"%s\" ;# Observation ID\n", s); } fprintf(obj_file, "#\n" "# List of objects used for cross-correlation\n" "# Input was: "); switch(cfg->saa_objacq_source) { case OBJACQ_AUTO: fprintf(obj_file, "automatic\n") ; break ; case OBJACQ_FILE: fprintf(obj_file, "user-provided file (%s)\n", cfg->saa_objacq_filename); break ; } fprintf(obj_file, "#\n" "#\tposx\tposy\n" "#\n" "\n" ); for (i=0 ; ixcorr_nobj ; i++) { fprintf(obj_file, "\t%d\t%d\n", cfg->xcorrobj_x[i], cfg->xcorrobj_y[i]); } fclose(obj_file); return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_find_offsets @memo Find out offsets between frames, relative to first frame. @param jit_cube Jitter cube. @param cfg Jitter blackboard. @return int 0 if Ok, anything else otherwise. @doc This function estimates the offsets in the cube, relative to the first frame. Offsets are stored in the blackboard. */ /*--------------------------------------------------------------------------*/ static int jitter_find_offsets(cube_t ** jit_cube, jitter_config * cfg) { double3 * offs, * estimates ; int i, j ; double3 * xcorr_obj ; cube_t * purged ; int * correct_frame ; int ncorrect ; double err_x, err_y ; if (cfg->saa_apply == 0) { e_comment(1, "no shift-and-add requested: skipping offset section") ; cfg->status_offacq = ALGO_SKIPPED ; return 0 ; } estimates = NULL ; switch (cfg->saa_offinput) { case OFFSIN_HEADER: /* Read header offsets */ estimates = get_offsets_from_fits_header(cfg) ; if (estimates==NULL) { /* Cannot get offsets from headers: exit now */ e_error("getting offsets from header"); cfg->status_offacq = ALGO_FAILED ; return -1 ; } break ; case OFFSIN_FILE: estimates = load_offsets_from_txtfile(cfg->saa_off_filename); if (estimates == NULL) { e_error("getting offsets from file"); cfg->status_offacq = ALGO_FAILED ; return -1 ; } else if (estimates->n < cfg->nobj) { e_error("not enough offset measurements found in file"); cfg->status_offacq = ALGO_FAILED ; return -1 ; } break ; case OFFSIN_AUTO: estimates = NULL ; break ; } /* Allocate space in the cfg structure to receive offsets */ cfg->xcorr_x = calloc(cfg->nobj, sizeof(double)) ; cfg->xcorr_y = calloc(cfg->nobj, sizeof(double)) ; cfg->xcorr_errx = calloc(cfg->nobj, sizeof(double)) ; cfg->xcorr_erry = calloc(cfg->nobj, sizeof(double)) ; cfg->xcorr_dist = calloc(cfg->nobj, sizeof(double)) ; /* * If some offsets have been received (estimates != NULL) * and refining has not been requested write the estimates into the cfg * struct, possibly dump offsets to file, and return. */ if ((estimates!=NULL) && (cfg->saa_off_refine==0)) { for (i=0 ; inobj ; i++) { cfg->xcorr_x[i] = estimates->x[i] ; cfg->xcorr_y[i] = estimates->y[i] ; } jitter_dump_offsets_file(cfg, estimates, NULL) ; double3_del(estimates); cfg->status_offacq = ALGO_OK ; return 0 ; } /* * If estimates have been received, check that they start at * (0,0), otherwise they are not in relative mode. If then, * subtract the first offsets from all consecutive offsets. */ if (estimates!=NULL) if ((X_ABS(estimates->x[0]) > 1) || (X_ABS(estimates->y[0]) > 1)) { e_warning("provided offsets do not start from (0,0) at 1st frame"); e_warning("subtracting first offsets in all frames"); for (i=1 ; inobj ; i++) { estimates->x[i] -= estimates->x[0] ; estimates->y[i] -= estimates->y[0] ; } estimates->x[0] = 0.0 ; estimates->y[0] = 0.0 ; } /* * Some work needs to be done, either because refining has been * requested, or because we are in automatic mode. * Test errors before launching cross-correlation. */ if ((cfg->saa_off_sx<0) || (cfg->saa_off_sy<0) || (cfg->saa_off_hx<0) || (cfg->saa_off_hy<0)) { e_error("in search or measure box"); e_error("search [%d %d] measure [%d %d]", cfg->saa_off_sx, cfg->saa_off_sy, cfg->saa_off_hx, cfg->saa_off_hy); cfg->status_offacq = ALGO_FAILED ; return -1 ; } /* * Retrieve stored info from cfg to know where xcorr objects are */ xcorr_obj = double3_new(cfg->xcorr_nobj); for (i=0 ; ixcorr_nobj ; i++) { xcorr_obj->x[i] = cfg->xcorrobj_x[i] ; xcorr_obj->y[i] = cfg->xcorrobj_y[i] ; xcorr_obj->z[i] = 0.0 ; } /* * Compute offsets between frames */ offs = xcorr_with_objs((*jit_cube), (*jit_cube)->plane[0], estimates, xcorr_obj, cfg->saa_off_sx, cfg->saa_off_sy, cfg->saa_off_hx, cfg->saa_off_hy); double3_del(xcorr_obj); if (offs==NULL) { e_error("estimating offsets") ; if (cfg->saa_off_refine==1 && estimates!=NULL) { e_warning("using estimates as offsets themselves") ; offs = estimates ; estimates = NULL ; } else { if (estimates!=NULL) double3_del(estimates) ; cfg->status_offacq = ALGO_FAILED ; return -1 ; } } /* Update blackboard */ cfg->xcorr_n = offs->n ; /* * Examine returned offsets to see if they make sense */ correct_frame = malloc(offs->n * sizeof(int)); ncorrect = 0 ; for (i=0 ; in ; i++) { if (estimates!=NULL) { err_x = fabs(offs->x[i] - estimates->x[i]) ; err_y = fabs(offs->y[i] - estimates->y[i]) ; } else { err_x = offs->x[i] ; err_y = offs->y[i] ; } err_x = fabs(err_x - (double)cfg->saa_off_sx); err_y = fabs(err_y - (double)cfg->saa_off_sy); if ((err_x < 0.1) || (err_y < 0.1) || (offs->z[i]<0.0)) { e_warning("object frame #%d cross-correlates badly", i+1); e_warning("removing this frame from the final stack"); correct_frame[i] = 0 ; offs->x[i] = 0.00 ; offs->y[i] = 0.00 ; offs->z[i] = -1.00 ; } else { correct_frame[i] = 1 ; ncorrect++ ; } } if (ncorrect<(offs->n/2)) { e_warning("less than half of the input frames correlate correctly"); switch (cfg->saa_offinput) { case OFFSIN_HEADER: e_warning("can we trust these header offset inputs?"); break ; case OFFSIN_FILE: e_warning("can we trust these file offset inputs?"); break ; case OFFSIN_AUTO: e_warning("the requested search size [%d x %d]", cfg->saa_off_sx, cfg->saa_off_sy); e_warning("is probably too small..."); break ; default: e_error("internal: in examining returned offsets"); if (offs!=NULL) double3_del(offs); if (estimates!=NULL) double3_del(estimates); free(correct_frame); return -1 ; } } if (ncorrect<1) { e_error("no frame correlated correctly: giving up"); if (offs!=NULL) double3_del(offs); if (estimates!=NULL) double3_del(estimates); free(correct_frame); return -1 ; } /* Update blackboard */ cfg->xcorr_nval = ncorrect ; /* Rearrange the cube if needed to drop mis-correlated planes */ if (ncorrectn) { purged = cube_new((*jit_cube)->lx, (*jit_cube)->ly, ncorrect); j=0 ; for (i=0 ; in ; i++) { if (correct_frame[i]) { purged->plane[j] = (*jit_cube)->plane[i]; j++ ; } else { image_del((*jit_cube)->plane[i]); (*jit_cube)->plane[i] = NULL ; } } cube_del_shallow(*jit_cube); *jit_cube = purged ; } free(correct_frame) ; cfg->nobj = ncorrect ; /* * Write offsets into jitter config */ jitter_dump_offsets_file(cfg, offs, estimates) ; double3_del(offs) ; if (estimates!=NULL) double3_del(estimates) ; cfg->status_offacq = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_shiftandadd @memo Apply the shift and add pass to a jitter cube. @param jit_cube Jitter cube. @param cfg Jitter blackboard @return int 0 if Ok, -1 otherwise. @doc Apply the shift-and-add process to a jitter cube. This modifies the input cube and returns it with a single image. */ /*--------------------------------------------------------------------------*/ static int jitter_shiftandadd(cube_t ** jit_cube, jitter_config * cfg) { int i, j ; double3 * offs ; image_t * final ; if (cfg->saa_apply!=1) return 0 ; /* Keep only correct offsets in list */ offs = double3_new(cfg->nobj); /* Copy valid offsets to offs */ j=0 ; for (i=0 ; i<(*jit_cube)->np; i++) { if (cfg->xcorr_dist[i]>-0.5) { offs->x[j] = - cfg->xcorr_x[i] ; offs->y[j] = - cfg->xcorr_y[i] ; offs->z[j] = cfg->xcorr_dist[i] ; j++ ; } } final = cube_shiftandadd( (*jit_cube), offs, NULL, cfg->saa_filt_lo, cfg->saa_filt_hi) ; double3_del(offs); if (final==NULL) { cfg->status_saa = ALGO_FAILED ; return -1 ; } cube_del(*jit_cube); (*jit_cube) = cube_new(final->lx, final->ly, 1); (*jit_cube)->plane[0] = final ; cfg->status_saa = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_postproc @memo Apply optional post-processings if required. @param jit_cube Jitter cube. @param cfg Jitter blackboard. @return int 0 if Ok, anything else otherwise. @doc This function applies optional post-processing to the jitter cube. */ /*--------------------------------------------------------------------------*/ static int jitter_postproc(cube_t * jit_cube, jitter_config * cfg) { int pn ; if (cfg->postproc_active==0) { e_comment(1, "post-processing deactivated: skipped") ; cfg->status_postproc = ALGO_SKIPPED ; return 0 ; } /* * median subtraction from each row */ if (cfg->postproc_rowmediansub == 0) { e_comment(1, "skipping: subtract the median from each row") ; } else { e_comment(1, "subtracting the median for each row") ; for (pn=0 ; pnnp ; pn++) { if (image_sub_rowmedian(jit_cube->plane[pn])!=0) { e_error("during median subtraction from each row") ; cfg->status_postproc = ALGO_FAILED ; return -1 ; } } } /* * If requested, output a QC report PAF file as basename_qc.paf */ if (cfg->postproc_qcreport) { if (jitter_qc(cfg, jit_cube) == -1) { e_error("cannot create QC PAF file") ; cfg->status_postproc = ALGO_FAILED ; return -1 ; } } /* * Other post-processing features here * ... */ cfg->status_postproc = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_qc @memo Quality control parameters derived from the output image @param cfg jitter configuration structure @param jit_cube jitter stacked output image (cube_t* type) @return 0 if ok, -1 otherwise @doc The created PAF file contains the following parameters: */ /*--------------------------------------------------------------------------*/ static int jitter_qc( jitter_config * cfg, cube_t * jit_cube) { FILE * paf ; char pafname[FILENAMESZ] ; char * strvar ; detected * det ; double fwhm_avg_med ; double pix_scale ; double iq ; int i ; e_comment(1, "creating output PAF file for QC1..."); sprintf(pafname, "%s_qc.paf", get_rootname(cfg->output_basename)) ; paf = paf_print_header( pafname, "ISAAC/jitter", "jitter recipe results") ; if (paf==NULL) { e_warning("cannot create PAF file: no QC output") ; return 0 ; } fprintf(paf, "\n"); /* MJD-OBS */ strvar = isaac_get_mjdobs(cfg->input_list[0]) ; if (strvar!=NULL) { fprintf(paf, "MJD-OBS %s; # Obs start\n\n", strvar); } else { fprintf(paf, "MJD-OBS 0.0; # Obs start unknown\n\n"); } /* INSTRUME keyword */ strvar = isaac_get_instrument(cfg->input_list[0]) ; if (strvar != NULL) { fprintf(paf, "INSTRUME \"%s\"\n", strvar) ; } /* TPL.ID */ strvar = isaac_get_templateid(cfg->input_list[0]) ; if (strvar != NULL) { fprintf(paf, "TPL.ID \"%s\"\n", strvar) ; } /* TPL.NEXP */ strvar = isaac_get_numbexp(cfg->input_list[0]) ; if (strvar != NULL) { fprintf(paf, "TPL.NEXP %s\n", strvar) ; } /* DPR.CATG */ strvar = isaac_get_dpr_catg(cfg->input_list[0]) ; if (strvar != NULL) { fprintf(paf, "DPR.CATG \"%s\"\n", strvar) ; } /* DPR.TYPE */ strvar = isaac_get_dpr_type(cfg->input_list[0]) ; if (strvar != NULL) { fprintf(paf, "DPR.TYPE \"%s\"\n", strvar) ; } /* DPR.TECH */ strvar = isaac_get_dpr_tech(cfg->input_list[0]) ; if (strvar != NULL) { fprintf(paf, "DPR.TECH \"%s\"\n", strvar) ; } /* Add PRO.CATG */ strvar = isaac_get_pro_catg_value(isaac_imag_sw_jitter_qc) ; if (strvar!=NULL) { fprintf(paf, "PRO.CATG \"%s\" ;# Product category\n", strvar); } /* Add the date */ strvar = isaac_get_date_obs(cfg->input_list[0]); if (strvar!=NULL) { fprintf(paf, "DATE-OBS \"%s\" ;# Date\n", strvar); } /* Add OBS.ID */ strvar = isaac_get_obs_id(cfg->input_list[0]); if (strvar!=NULL) { fprintf(paf, "OBS.ID %s ;# Obs id\n", strvar); } /* Add INS.PIXSCALE */ strvar = isaac_get_pixscale(cfg->input_list[0]) ; if (strvar != NULL) { pix_scale = (double)atof(strvar); fprintf(paf, "INS.PIXSCALE %s\n", strvar) ; } else { pix_scale = -1.0 ; } /* Write the stats in the PAF file */ /* Write out sky background measurements if any */ if (cfg->sky_background!=NULL) { e_comment(2, "printing out sky background measurements..."); fprintf(paf, "\n\nJITTER.SKYBG.START\n"); for (i=0 ; inframes ; i++) { fprintf(paf, "%g\n", cfg->sky_background[i]); } fprintf(paf, "JITTER.SKYBG.END\n\n\n"); } /* Compute statistics on output jittered frame */ e_comment(2, "detecting objects on final frame..."); if ((det = detected_ks_engine(jit_cube->plane[0], DETECTED_KAPPA, 0)) == NULL) { e_warning("cannot find objects on result frame...") ; fwhm_avg_med = -1.0 ; } else { e_comment(2, "computing median FWHM on final frame..."); /* Compute FWHMs */ detected_compute_fwhm(det, jit_cube->plane[0]); /* Display results in PAF file */ fprintf(paf, "JITTER.OBJECTS.START\n") ; detected_dump(det, paf); fprintf(paf, "JITTER.OBJECTS.END\n") ; /* Save median FWHM in frame */ fwhm_avg_med = det->fwhm_meda ; /* Compute image quality */ iq = detected_compute_iq(det, pix_scale, NULL); detected_del(det) ; } fprintf(paf, "\n\n"); if (fwhm_avg_med>0.0) { e_comment(2, "median FWHM: %g pixels", fwhm_avg_med); fprintf(paf, "QC.FWHM.PIX %g \n", fwhm_avg_med) ; } /* Median FWHM in arcsec */ if (pix_scale>0.0) { e_comment(2, "median FWHM: %g arcsec", fwhm_avg_med * pix_scale); fprintf(paf, "QC.FWHM.ARCSEC %g\n", fwhm_avg_med * pix_scale) ; } /* Image quality in arcsec */ if (iq>0.0) { e_comment(2, "image quality: %g arcsec", iq); fprintf(paf, "QC.IQ %g\n", iq); } /* Add FILTER */ strvar = isaac_get_filter(cfg->input_list[0]); if (strvar!=NULL) { fprintf(paf, "QC.FILTER.OBS \"%s\"\n", strvar); } e_comment(1, "output PAF file complete"); fclose(paf); return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_save_cube @memo Save the output of the jitter process. @param jit_cube Jitter cube. @param cfg Jitter blackboard. @return int 0 if Ok, anything else otherwise. @doc Save the output of the jitter process to a FITS file. */ /*--------------------------------------------------------------------------*/ static int jitter_save_cube(cube_t * jit_cube, jitter_config * cfg) { char outname[1024] ; qfits_header* fh ; framelist * input_files ; /* Read FITS header from the first frame */ fh = qfits_header_read(cfg->input_list[0]); if (fh==NULL) { e_error("cannot read header from %s: creating empty header", cfg->input_list[0]); } /* Make header fit for an image */ if (fh!=NULL) { if (isaac_header_for_image(fh)!=0) { e_error("cleaning input header: creating empty header"); qfits_header_destroy(fh); fh=NULL ; } } /* Update FITS header with PRO keywords */ sprintf(outname, "%s.fits", cfg->output_basename); input_files = framelist_load(cfg->frames_filelist) ; isaac_pro_fits(fh, outname, "REDUCED", "PHOTOMETRIC", isaac_imag_sw_jitter_result, "OK", "img_jitter", cfg->nobj, input_files, NULL); framelist_del(input_files) ; /* Add various comments */ jitter_add_comments(cfg, fh) ; e_comment(0, "saving final output [%s]", outname); cube_save_fits_hdrdump(jit_cube, outname, fh); qfits_header_destroy(fh); cfg->status_save = ALGO_OK ; return 0 ; } /*--------------------------------------------------------------------------- Function : jitter_offset_skysub() In : cube to process, SYMTAB Out : int 0 of Ok, anything else otherwise Job : remove the sky from a batch of jitter+offset frames Notice : poor algorithm... ---------------------------------------------------------------------------*/ static int jitter_offset_skysub(cube_t ** jit, jitter_config * cfg) { int nframes ; int n_obj, n_sky ; image_t ** all_obj ; cube_t * all_sky ; cube_t * skysub ; image_t * median_sky ; pixelvalue median_pix ; char outputdiff_name[FILENAMESZ] ; int nbpix ; qfits_header * hdr ; int i, j, k ; /* Count how many frames of each type there are */ nframes = cfg->nframes ; n_obj = cfg->nobj ; n_sky = nframes - n_obj ; /* * Store all frames into two structures */ all_obj = malloc(n_obj * sizeof (image_t*)); all_sky = cube_new((*jit)->lx, (*jit)->ly, n_sky); j=0 ; k=0 ; for (i=0 ; iframetype[i]==FRAME_SKY) { all_sky->plane[j++] = (*jit)->plane[i] ; } else { all_obj[k++] = (*jit)->plane[i] ; } } e_comment(2, "computing median sky...") ; median_sky = cube_avg_median(all_sky) ; free(all_sky->plane); free(all_sky); /* * Now subtract median sky from all input object planes */ e_comment(2, "subtracting median sky from all objects..."); skysub = cube_new((*jit)->lx, (*jit)->ly, n_obj); for (i=0 ; iplane[i] = image_sub(all_obj[i], median_sky); } image_del(median_sky); free(all_obj); cube_del(*jit); /* * Now subtract median value from each object plane */ nbpix = skysub->plane[0]->lx * skysub->plane[0]->ly ; for (i=0 ; iplane[i]!=NULL) { median_pix = image_getmedian(skysub->plane[i]); for (j=0 ; jplane[i]->data[j] -= median_pix ; } } } if (cfg->sky_outputdiff) { e_comment(1, "saving sky-subtracted frames"); /* Read the fits header of the first input frame */ if ((hdr = qfits_header_read(cfg->input_list[0])) == NULL) { e_error("cannot output sky bg frame") ; } else { sprintf(outputdiff_name, "%s_dif.fits", cfg->output_basename); isaac_header_for_image(hdr) ; isaac_pro_fits(hdr, outputdiff_name, "REDUCED", NULL, isaac_imag_sw_jitter_diff, "OK", "img_jitter", cfg->nframes, NULL, NULL) ; cube_save_fits_hdrdump(skysub, outputdiff_name, hdr) ; qfits_header_destroy(hdr) ; } } (*jit) = skysub ; cfg->status_skyengine = ALGO_OK ; cfg->status_skybg = ALGO_SKIPPED ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_add_comments @memo Add comments to the jitter output FITS file. @param cfg Jitter blackboard. @param fh FITS header to update. @return void @doc This function fetches various processing informations from the blackboard and adds them as commens to the output jitter cube. These will end up as HISTORY keywords in the output FITS file. */ /*--------------------------------------------------------------------------*/ static void jitter_add_comments(jitter_config * cfg, qfits_header * fh) { history * hs ; int i ; /* * The following data will be written in the output header as * HISTORY fields: * * Version of the jitter software * List of input files * List of flat-field and dark frames if used * Sky estimation statistics (background measurements) * Which object was used for Xcorrelation * Offset measurements + reliability * Total exposure time */ hs = history_new(); history_add(hs, "--- eclipse jitter imaging data reduction"); history_add(hs, "jitter software version: %s", get_eclipse_version()) ; history_add(hs, "[AlgorithmStatus]"); history_add(hs, "Cosmetics = %s", jitter_algo(cfg->status_cosmetics)); history_add(hs, "SkyEngine = %s", jitter_algo(cfg->status_skyengine)); history_add(hs, "SkyBackGround = %s", jitter_algo(cfg->status_skybg)); history_add(hs, "ObjectAcq = %s", jitter_algo(cfg->status_objacq)); history_add(hs, "OffsetAcq = %s", jitter_algo(cfg->status_offacq)); history_add(hs, "ShiftAndAdd = %s", jitter_algo(cfg->status_saa)); history_add(hs, "PostProc = %s", jitter_algo(cfg->status_postproc)); history_add(hs, "[Frames]"); history_add(hs, "FileList = %s", cfg->frames_filelist); history_add(hs, "Path = %s", get_dirname(cfg->input_list[0])) ; history_add(hs, "NFrames = %d", cfg->nframes); history_add(hs, "NObj = %d", cfg->nobj) ; for (i=0 ; inframes ; i++) { history_add(hs, "Frame:%03d (%s) %s", i+1, IFTYP(cfg->frametype[i]), get_basename(cfg->input_list[i])); } history_add(hs, "[Dark]"); history_add(hs, "Subtraction = %s", IBOOL(cfg->dark_subtraction)); if (cfg->dark_subtraction) { history_add(hs, "Filename = %s", get_basename(cfg->dark_filename)); } history_add(hs, "[FlatField]"); history_add(hs, "Division = %s", IBOOL(cfg->flatfield_division)); if (cfg->flatfield_division) { history_add(hs, "Filename = %s", get_basename(cfg->flatfield_filename)); } history_add(hs, "[BadPixels]"); history_add(hs, "Replacement = %s", IBOOL(cfg->badpixels_replacement)); if (cfg->badpixels_replacement) { history_add(hs, "Filename = %s", get_basename(cfg->badpixels_filename)); } history_add(hs, "[SkyEngine]"); history_add(hs, "EstimateSky = %s", IBOOL(cfg->sky_estimate)); if (cfg->sky_estimate) { history_add(hs, "MinNumberOfFrames = %d", cfg->sky_minnumofframes); history_add(hs, "[SkyCombine]"); history_add(hs, "RejectHalfWidth = %d", cfg->sky_rejhalfwidth); history_add(hs, "RejectMin = %d", cfg->sky_rejmin); history_add(hs, "RejectMax = %d", cfg->sky_rejmax); } history_add(hs, "[ShiftAndAdd]"); history_add(hs, "ApplyShiftAndAdd = %s", IBOOL(cfg->saa_apply)); history_add(hs, "[PostProcessing]"); if (cfg->postproc_active && cfg->postproc_rowmediansub) { history_add(hs, "RowSubtractMedian = yes"); } else { history_add(hs, "RowSubtractMedian = no"); } /* Dump history into FITS header */ history_addfits(hs, fh); history_del(hs); return ; }