/*--------------------------------------------------------------------------- File name : lwjit_engine.c Author : N. Devillard Created on : July 2000 Description : lwjitter main engine ---------------------------------------------------------------------------*/ /* $Id: lwjit_engine.c,v 1.33 2001/10/22 11:57:22 ndevilla Exp $ $Author: ndevilla $ $Date: 2001/10/22 11:57:22 $ $Revision: 1.33 $ */ /*--------------------------------------------------------------------------- Includes ---------------------------------------------------------------------------*/ #include #include #include #include #include #include "eclipse.h" #include "lwjit_bb.h" #include "lwjit_engine.h" #include "lwjit_offset.h" #include "lwjit_class.h" #include "lwjit_ini.h" #include "isaacp_lib.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") /*--------------------------------------------------------------------------- Functions private to this module ---------------------------------------------------------------------------*/ static cube_t* lwjit_subtract_chopnod(lwjit_bb*); static int jitter_acquire_objects(cube_t *, lwjit_bb *); static int jitter_get_object_source(lwjit_bb *); static int jitter_load_objects_from_file(lwjit_bb *); static int jitter_autofindobjects(cube_t *, lwjit_bb *); static int jitter_output_object_list(lwjit_bb *); static int jitter_find_offsets(cube_t **, lwjit_bb *); static void switch_auto_find_offsets(lwjit_bb *); static int jitter_register_planes(cube_t *, lwjit_bb *) ; static int jitter_stack_frames(cube_t **, lwjit_bb *); static int jitter_postproc(cube_t *, lwjit_bb *) ; static int jitter_save_cube(cube_t * jit_cube, lwjit_bb *); /*---------------------------------------------------------------------------*/ /** @name lwjit_engine @memo Launch the lwjitter engine from a given configuration. @param ini_name Name of ini file to use to initialize the engine @return long number of pixels processed in input. @doc Main lwjitter engine function. The whole algorithm can be parsed here. */ /*---------------------------------------------------------------------------*/ #define ALGO_PARTS 9 long lwjit_engine(char * ini_name) { lwjit_bb * bb ; cube_t * jit_cube ; time_t local_t ; char filename[FILENAMESZ] ; char status_file[FILENAMESZ] ; long total_pixels_in ; int err ; int part ; part=0 ; e_comment(0, "parsing ini file [%s]", ini_name); bb = lwjit_ini_parse(ini_name); if (bb==NULL) { printf("error parsing file [%s]: aborting\n", ini_name); return -1 ; } e_comment(0, "--> START lwjit engine\n"); time(&local_t) ; e_comment(0, "%s", ctime(&local_t)) ; e_comment(0, "pid is %ld\n", (long)getpid()); /* Classification of input frames */ part++ ; e_comment(0, "---> part %d of %d: classification of frames", part, ALGO_PARTS) ; err = lwjit_classification(bb) ; if (err!=0) { e_error("cannot classify input frames") ; lwjit_bb_destroy(bb); return -1 ; } part++; e_comment(0, "---> part %d of %d: chop/nod subtraction", part, ALGO_PARTS) ; jit_cube=lwjit_subtract_chopnod( bb ); if (jit_cube==NULL) { e_error("in chopnod subtraction: skipping") ; lwjit_bb_destroy(bb); return -1; } total_pixels_in = 2 * (long)(jit_cube->lx * jit_cube->ly * jit_cube->np) ; /* Apply dark subtraction/flat-field division if requested */ part++; e_comment(0, "---> part %d of %d: flat-field & bad pixels", part, ALGO_PARTS) ; if ((bb->flatfield_division == 0) && (bb->badpixels_replacement == 0)) { bb->status_cosmetics = ALGO_SKIPPED ; } else { isaac_ff_dark_badpix_handling(&jit_cube, bb->flatfield_filename, NULL, bb->badpixels_filename); bb->status_cosmetics = ALGO_OK ; } /* Find out objects in first frame, for x-correlation */ part++; e_comment(0, "---> part %d of %d: object acquisition", part, ALGO_PARTS); if (jitter_acquire_objects(jit_cube, bb)!=0) { e_error("in object acquisition: aborting"); cube_del(jit_cube); lwjit_bb_destroy(bb); return -1 ; } /* Find out the offsets */ part++; e_comment(0, "---> part %d of %d: offset determination", part, ALGO_PARTS) ; if (jitter_find_offsets(&jit_cube, bb)!=0) { e_error("in offset finding section: aborting") ; cube_del(jit_cube) ; lwjit_bb_destroy(bb); return -1 ; } /* Apply the offsets to all planes */ part++; e_comment(0, "---> part %d of %d: plane registration", part, ALGO_PARTS) ; if (jitter_register_planes(jit_cube, bb)!=0) { e_error("in plane registration") ; cube_del(jit_cube); lwjit_bb_destroy(bb); return -1 ; } /* Average the result to a single plane */ part++; e_comment(0, "---> part %d of %d: 3d average", part, ALGO_PARTS); if (jitter_stack_frames(&jit_cube, bb)!=0) { e_error("in frame stacking") ; } /* Optional post-processing features will be inserted here */ part++; e_comment(0, "---> part %d of %d: optional post-processing", part, ALGO_PARTS) ; if (jitter_postproc(jit_cube, bb)!=0) { e_error("in post-processing") ; } /* Save results */ part++; e_comment(0, "---> part %d of %d: saving output data", part, ALGO_PARTS); jitter_save_cube(jit_cube, bb) ; cube_del(jit_cube) ; e_comment(0, "--> STOP lwjit engine\n"); time(&local_t) ; e_comment(0, "%s", ctime(&local_t)) ; /* If requested: launch an image viewer on the result */ if ((bb->postproc_startviewer == 1) && (bb->postproc_active == 1)) { sprintf(filename, "%s.fits", bb->output_basename); show_image(filename, bb->postproc_viewer) ; } /* If requested, output a status file as basename_status.ascii */ if (bb->postproc_statusreport) { sprintf(status_file, "%s_status.ascii", bb->output_basename); lwjit_bb_dump(bb, status_file); } lwjit_bb_destroy(bb); return total_pixels_in ; } /*-------------------------------------------------------------------------*/ /** @name lwjit_subtract_chopnod @memo Estimate and subtract the nod frames. @param bb Jitter blackboard. @return int 0 if Ok, -1 otherwise. @doc Estimates and subtract the nod frames. */ /*--------------------------------------------------------------------------*/ static cube_t *lwjit_subtract_chopnod( lwjit_bb *bb) { cube_t * chop_a ; cube_t * chop_b ; chop_a = cube_load_strings(bb->chop_a, bb->nnod); if (chop_a == NULL) { e_error("loading first batch of frames: aborting"); bb->status_subtraction = ALGO_FAILED ; return NULL ; } chop_b = cube_load_strings(bb->chop_b, bb->nnod); if (chop_b == NULL) { e_error("loading second batch of frames: aborting"); cube_del(chop_a); bb->status_subtraction = ALGO_FAILED ; return NULL ; } e_comment(1, "subtracting pairs..."); cube_sub(chop_a, chop_b); cube_del(chop_b); cube_cst_op(chop_a, 0.5, '*'); bb->status_subtraction = ALGO_OK ; return chop_a; } /*-------------------------------------------------------------------------*/ /** @name jitter_acquire_objects @memo Find object for x-correlation. @param jit_cube Jitter cube @param bb Jitter blackboard @return int 0 if Ok, anything else otherwise. @doc Acquire a list of objects for cross-correlation. Objects are either acquired from an ASCII file provided by the user, or automatically through an object detection algorithm. */ /*--------------------------------------------------------------------------*/ static int jitter_acquire_objects(cube_t * jit_cube, lwjit_bb * bb) { int obj_source ; int status ; obj_source = jitter_get_object_source(bb); switch (obj_source) { case XCORR_SRC_ERROR: status = -1 ; bb->status_objacq = ALGO_FAILED ; break ; case XCORR_SRC_NOACQ: e_comment(1, "no need to acquire objects for x-correlation"); bb->status_objacq = ALGO_SKIPPED ; return 0 ; case XCORR_SRC_AUTO: e_comment(1, "automatic object acquisition for x-correlation"); status = jitter_autofindobjects(jit_cube, bb); break ; case XCORR_SRC_FILE: e_comment(1, "loading objects from user source"); status = jitter_load_objects_from_file(bb); break; default: status = -1 ; break ; } if (status==0) { bb->status_objacq = ALGO_OK ; jitter_output_object_list(bb); } else bb->status_objacq = ALGO_FAILED ; return status ; } /*-------------------------------------------------------------------------*/ /** @name jitter_get_object_source @memo Find where x-correlation objects are coming from. @param bb Jitter blackboard. @return int code, see below. @doc This function looks up the jitter blackboard to determine where objects are expected from. It returns one of the following codes: \begin{tabular}{ll} XCORR_SRC_ERROR & An error occurred \\ XCORR_SRC_NOACQ & No objects need to be acquired \\ XCORR_SRC_AUTO & Acquisition in automatic mode \\ XCORR_SRC_FILE & Acquisition done from file \end{tabular} */ /*--------------------------------------------------------------------------*/ static int jitter_get_object_source(lwjit_bb * bb) { int need_acq ; int need_detect ; /* * 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. */ if (bb->saa_apply==0) { e_comment(1, "no shift and add requested: skipping object acq"); return XCORR_SRC_NOACQ ; } need_acq = 0 ; switch (bb->off_input) { case OFFSIN_HEADER: need_acq = bb->off_hdr_refine ; break ; case OFFSIN_FILE: need_acq = bb->off_file_refine ; break ; case OFFSIN_AUTO: need_acq = 1 ; break ; default: e_error("internal: illegal offset input type: %d", bb->off_input); return XCORR_SRC_ERROR ; } if (need_acq==0) { e_comment(1, "no refining requested: no need to detect objects"); return XCORR_SRC_NOACQ ; } /* * Now it is sure that objects must be acquired. Find out where the * source comes from. */ switch (bb->objacq_source) { case OBJACQ_AUTO: need_detect = 1 ; break ; case OBJACQ_FILE: need_detect = 0 ; break ; default: e_warning("internal: illegal object source: [%d]",bb->objacq_source); e_warning("switching to automatic object finder"); need_detect = 0 ; } if (need_detect == 0) return XCORR_SRC_FILE ; return XCORR_SRC_AUTO ; } /*-------------------------------------------------------------------------*/ /** @name jitter_autofindobjects @memo Find automatically relevant objects for x-correlation. @param jit_cube Jitter cube. @param bb 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, lwjit_bb * bb) { double3 * peaks ; int edge_x, edge_y ; int i ; /* Get object detection parameters */ switch (bb->off_input) { case OFFSIN_HEADER: edge_x = bb->off_hdr_sx + bb->off_hdr_hx ; edge_y = bb->off_hdr_sy + bb->off_hdr_hy ; break ; case OFFSIN_FILE: edge_x = bb->off_file_sx + bb->off_file_hx ; edge_y = bb->off_file_sy + bb->off_file_hy ; break ; case OFFSIN_AUTO: edge_x = bb->off_auto_sx + bb->off_auto_hx ; edge_y = bb->off_auto_sy + bb->off_auto_hy ; break ; default: e_error("unknown detection parameter off_input=%d",bb->off_input); return -1; } /* Launch object detection parameters */ peaks = get_xcorrelation_points(jit_cube->plane[0], edge_x, edge_y, bb->detpeak_threshold, bb->detpeak_mindetpoints, bb->detpeak_maxdetpoints); if (peaks==NULL) { e_error("cannot find enough valid points for xcorrelation"); return -1 ; } /* Update relevant fields in config table */ bb->xcorr_nobj = peaks->n ; bb->xcorrobj_x = malloc(peaks->n * sizeof(int)) ; bb->xcorrobj_y = malloc(peaks->n * sizeof(int)) ; for (i=0 ; in ; i++) { bb->xcorrobj_x[i] = (int)peaks->x[i] ; bb->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 bb 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(lwjit_bb *bb) { char * obj_list_name ; double3 * obj_list ; int i ; /* Get the name of the object list from the bb table */ obj_list_name = bb->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 */ bb->xcorr_nobj = obj_list->n ; bb->xcorrobj_x = malloc(obj_list->n*sizeof(int)); bb->xcorrobj_y = malloc(obj_list->n*sizeof(int)); for (i=0 ; in ; i++) { bb->xcorrobj_x[i] = (int)(obj_list->x[i]+0.5) ; bb->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 bb 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(lwjit_bb * bb) { char obj_file_name[FILENAMESZ] ; FILE * obj_file ; int i ; if (bb->detpeak_outputobjs == 0) return 0 ; sprintf(obj_file_name, "%s_obj.ascii", bb->output_basename); if ((obj_file=fopen(obj_file_name, "w"))==NULL) { e_error("cannot create [%s]: aborting output", obj_file_name); return -1 ; } if (bb->xcorr_nobj<1) { e_error("strange number of objects [%d]: aborting output", bb->xcorr_nobj); fclose(obj_file); return -1 ; } fprintf(obj_file, "#\n# List of objects used for cross-correlation\n# Input was: "); switch(bb->objacq_source) { case OBJACQ_AUTO: fprintf(obj_file, "automatic\n") ; break ; case OBJACQ_FILE: fprintf(obj_file, "user-provided file (%s)\n", bb->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", bb->xcorrobj_x[i], bb->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 bb 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, lwjit_bb * bb) { double3 * offs, * estimates ; int i, j ; int refine_flag ; int sx, sy ; int hx, hy ; double3 * xcorr_obj ; cube_t * purged ; int * correct_frame ; int ncorrect ; double err_x, err_y ; if (bb->saa_apply == 0) { e_comment(1, "no shift-and-add requested: skipping offset section") ; bb->status_offacq = ALGO_SKIPPED ; return 0 ; } refine_flag = 0 ; estimates = NULL ; switch (bb->off_input) { case OFFSIN_HEADER: /* read header offsets */ estimates = get_offsets_from_fits_header(bb) ; if (estimates==NULL){ e_warning("cannot get offsets from header: switching to auto"); switch_auto_find_offsets(bb) ; refine_flag = 0 ; } else { refine_flag = bb->off_hdr_refine ; } break ; case OFFSIN_FILE: estimates = load_offsets_from_txtfile(bb->off_file_name); if (estimates == NULL) { e_warning("cannot get offsets from file: switching to auto"); estimates = NULL ; switch_auto_find_offsets(bb) ; refine_flag = 0 ; } else if (estimates->n < bb->nnod) { e_warning("not enough offsets in file: switching to auto"); double3_del(estimates); estimates = NULL ; switch_auto_find_offsets(bb); refine_flag = 0 ; } else { refine_flag = bb->off_file_refine ; } break ; case OFFSIN_AUTO: refine_flag = 0 ; estimates = NULL ; break ; } /* Allocate space in the bb structure to receive offsets */ bb->xcorr_x = calloc(bb->nnod, sizeof(double)) ; bb->xcorr_y = calloc(bb->nnod, sizeof(double)) ; bb->xcorr_errx = calloc(bb->nnod, sizeof(double)) ; bb->xcorr_erry = calloc(bb->nnod, sizeof(double)) ; bb->xcorr_dist = calloc(bb->nnod, sizeof(double)) ; /* * If we received some offsets (estimates != NULL) and refining has not * been requested (refine_flag == FALSE), write the estimates into the * bb struct, possibly dump offsets to file, and return. */ if ((estimates!=NULL) && (refine_flag==0)) { for (i=0 ; innod ; i++) { bb->xcorr_x[i] = estimates->x[i] ; bb->xcorr_y[i] = estimates->y[i] ; } jitter_dump_offsets_file(bb, estimates, NULL) ; double3_del(estimates); bb->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 ; innod ; 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. */ switch (bb->off_input) { case OFFSIN_HEADER: sx = bb->off_hdr_sx ; sy = bb->off_hdr_sy ; hx = bb->off_hdr_hx ; hy = bb->off_hdr_hy ; break ; case OFFSIN_FILE: sx = bb->off_file_sx ; sy = bb->off_file_sy ; hx = bb->off_file_hx ; hy = bb->off_file_hy ; break ; default: sx = bb->off_auto_sx ; sy = bb->off_auto_sy ; hx = bb->off_auto_hx ; hy = bb->off_auto_hy ; break ; } /* Test errors before launching cross-correlation */ if ((sx<0) || (sy<0) || (hx<0) || (hy<0)) { e_error("in search or measure box"); e_error("search [%d %d] measure [%d %d]", sx, sy, hx, hy); bb->status_offacq = ALGO_FAILED ; return -1 ; } /* Retrieve stored info from bb to know where xcorr objects are */ xcorr_obj = double3_new(bb->xcorr_nobj); for (i=0 ; ixcorr_nobj ; i++) { xcorr_obj->x[i] = bb->xcorrobj_x[i] ; xcorr_obj->y[i] = bb->xcorrobj_y[i] ; xcorr_obj->z[i] = (pixelvalue)0; } /* * Compute offsets between frames */ offs = xcorr_with_objs((*jit_cube), (*jit_cube)->plane[0], estimates, xcorr_obj, sx, sy, hx, hy); double3_del(xcorr_obj); if (offs==NULL) { e_error("getting offsets") ; if (refine_flag==1 && estimates!=NULL) { e_warning("using estimates as offsets themselves") ; offs = estimates ; estimates = NULL ; } else { if (estimates!=NULL) double3_del(estimates) ; bb->status_offacq = ALGO_FAILED ; return -1 ; } } /* Update blackboard */ bb->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)sx); err_y = fabs(err_y - (double)sx); 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 (bb->off_input) { 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]", sx, 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 */ bb->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) ; bb->nnod = ncorrect ; /* * Write offsets into jitter config */ jitter_dump_offsets_file(bb, offs, estimates) ; double3_del(offs) ; if (estimates!=NULL) double3_del(estimates) ; bb->status_offacq = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name switch_auto_find_offsets @memo Switch blackboard config to automatic offset finding. @param bb Jitter blackboard. @return void @doc This function acts upon the jitter blackboard to switch its configuration to automatic offset finding. This is needed when e.g. the provided offsets are incorrect. */ /*--------------------------------------------------------------------------*/ static void switch_auto_find_offsets(lwjit_bb * bb) { bb->off_input = OFFSIN_AUTO ; /* * If no value was provided for [sx,sy] and [hx,hy], * provide default ones. */ if (bb->off_auto_sx<1) bb->off_auto_sx = OFFS_AUTO_SX ; if (bb->off_auto_sy<1) bb->off_auto_sx = OFFS_AUTO_SY ; if (bb->off_auto_hx<1) bb->off_auto_sx = OFFS_AUTO_MX ; if (bb->off_auto_hy<1) bb->off_auto_sx = OFFS_AUTO_MY ; return ; } /*-------------------------------------------------------------------------*/ /** @name jitter_register_planes @memo Register all planes in the cube according to the offsets. @param jit_cube Jitter cube. @param bb Jitter blackboard. @return int 0 if Ok, anything else otherwise. @doc This function modifies the jitter cube to register all planes to the same reference: the first plane. Offsets are taken from the blackboard. */ /*--------------------------------------------------------------------------*/ static int jitter_register_planes(cube_t * jit_cube, lwjit_bb * bb) { int i, j ; double3 * offs ; int ret ; if (bb->saa_apply==0) { e_comment(1, "no shift-and-add requested: skipping register section") ; bb->status_register = ALGO_SKIPPED ; return 0 ; } /* * Get valid offsets into a double3 */ offs = double3_new(bb->xcorr_nval); j=0 ; for (i=0 ; ixcorr_nval ; i++) { if (bb->xcorr_dist[i]>-0.5) { offs->x[j] = bb->xcorr_x[i] ; offs->y[j] = bb->xcorr_y[i]; offs->z[j] = bb->xcorr_dist[i] ; j++ ; } } /* * Plane registration */ if (bb->reg_subpix) { ret = cube_shift(jit_cube, offs, bb->reg_kernel); } else { ret = cube_shift_int(jit_cube, offs) ; } double3_del(offs) ; if (jit_cube->np<1) { e_error("no valid plane for cross-correlation found in cube"); bb->status_register = ALGO_FAILED ; return 1 ; } if (ret!=0) { e_error("during cube shift: aborting shift") ; bb->status_register = ALGO_FAILED ; return 1 ; } bb->status_register = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_stack_frames @memo Stack frames in the jitter cube to a single one. @param jit_cube Jitter cube. @param bb Jitter blackboard. @return int 0 if Ok, 1 otherwise. @doc This function stacks the frames in the jitter cube to a single frame. It does nothing if stacking has not been requested. */ /*--------------------------------------------------------------------------*/ static int jitter_stack_frames(cube_t ** jit_cube, lwjit_bb * bb) { image_t * avg ; /* Average the result if required */ if (bb->saa_apply == 0) { e_comment(1, "no shift and add: skipped averaging") ; bb->status_3davg = ALGO_SKIPPED ; return 0 ; } if (bb->avg3d_active == 0) { e_comment(1, "skipped averaging") ; bb->status_3davg = ALGO_SKIPPED ; return 0 ; } if ((bb->avg3d_type == AVG3D_MEDIAN) || (bb->avg3d_type == AVG3D_FILTERED)) { if ((*jit_cube)->np<3) { e_warning("only [%d] frames in output stack", (*jit_cube)->np); e_warning("cannot do median or filtered 3d average"); e_warning("performing linear median"); bb->avg3d_type = AVG3D_LINEAR ; } } switch (bb->avg3d_type) { case AVG3D_LINEAR: avg = cube_avg_linear(*jit_cube) ; break ; case AVG3D_MEDIAN: avg = cube_avg_median(*jit_cube) ; break ; case AVG3D_FILTERED: avg = cube_avg_reject(*jit_cube, (int)(bb->avg3d_filt_lo * (*jit_cube)->np), (int)(bb->avg3d_filt_hi * (*jit_cube)->np)); break ; default: e_error("unrecognized average type: skipping average") ; bb->status_3davg = ALGO_SKIPPED ; return 1 ; } if (avg==NULL) { e_error("error in averaging cube: output will be a cube") ; bb->status_3davg = ALGO_FAILED ; return 1 ; } cube_del(*jit_cube) ; *jit_cube = cube_from_image(avg) ; image_del(avg) ; bb->status_3davg = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_postproc @memo Apply optional post-processings if required. @param jit_cube Jitter cube. @param bb 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, lwjit_bb * bb) { if (bb->postproc_active==0) { e_comment(1, "post-processing deactivated: skipped") ; bb->status_postproc = ALGO_SKIPPED ; return 0 ; } /* * No post-processing implemented yet */ bb->status_postproc = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name jitter_save_cube @memo Save the output of the jitter process. @param jit_cube Jitter cube. @param bb 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, lwjit_bb * bb) { char outname[1024] ; qfits_header* fh ; framelist * rawfiles ; /* Read FITS header from the first frame */ fh = qfits_header_read(bb->chop_a[0]); isaac_header_for_image(fh) ; /* Update FITS header with PRO keywords */ sprintf(outname, "%s.fits", bb->output_basename); rawfiles = framelist_load(bb->frames_filelist) ; isaac_pro_fits(fh, outname, "REDUCED", "PHOTOMETRIC", isaac_imag_lw_jitter_result, "OK", "isaaclw_img_jitter", bb->nnod, rawfiles, NULL); framelist_del(rawfiles) ; e_comment(0, "saving final output [%s]", outname); cube_save_fits_hdrdump(jit_cube, outname, fh); qfits_header_destroy(fh); bb->status_save = ALGO_OK ; return 0 ; }