sinfo_utl_illumcorr.c

00001 /* $Id: sinfo_utl_illumcorr.c,v 1.11 2006/12/19 17:01:36 amodigli Exp $
00002  *
00003  * This file is part of the CPL (Common Pipeline Library)
00004  * Copyright (C) 2002 European Southern Observatory
00005  *
00006  * This library is free software; you can redistribute it and/or
00007  * modify it under the terms of the GNU Lesser General Public
00008  * License as published by the Free Software Foundation; either
00009  * version 2.1 of the License, or (at your option) any later version.
00010  *
00011  * This library is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with this library; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 /*
00021  * $Author: amodigli $
00022  * $Date: 2006/12/19 17:01:36 $
00023  * $Revision: 1.11 $
00024  * $Name:  $
00025  */
00026 
00027  /****************************************************************
00028   *           Object Data reduction                              *
00029   ****************************************************************/
00030 
00031 #ifdef HAVE_CONFIG_H
00032 #include <config.h>          /* allows the program compilation */
00033 #endif
00034 
00035 /*-----------------------------------------------------------------------------
00036                                 Includes
00037 -----------------------------------------------------------------------------*/
00038 
00039 /* std */
00040 #include <strings.h>
00041 #include <string.h>
00042 #include <stdio.h>
00043 #include <math.h>
00044 #include <libgen.h>
00045 
00046 
00047 /* cpl */
00048 #include <cpl.h>  
00049 
00050 /* irplib */
00051 #include <irplib_utils.h>
00052 
00053 /* sinfoni */
00054 #include <sinfo_pro_types.h>
00055 #include <sinfo_product_config.h>
00056 #include <sinfo_prepare_stacked_frames_config.h>
00057 #include <sinfo_objnod_config.h>
00058 #include <sinfo_new_objnod.h>
00059 #include <sinfo_new_prepare_stacked_frames.h>
00060 #include <sinfo_pro_save.h>
00061 #include <sinfo_raw_types.h>
00062 #include <sinfo_functions.h>
00063 #include <sinfo_tpl_utils.h>
00064 #include <sinfo_tpl_dfs.h>
00065 #include <sinfo_hidden.h>
00066 #include <sinfo_globals.h>
00067 #include <sinfo_rec_utils.h>
00068 
00069 
00070 /*-----------------------------------------------------------------------------
00071                             Function prototypes
00072 -----------------------------------------------------------------------------*/
00073 static int sinfo_utl_illumcorr_create(cpl_plugin *plugin);
00074 static int sinfo_utl_illumcorr_exec(cpl_plugin *plugin);
00075 static int sinfo_utl_illumcorr_destroy(cpl_plugin *plugin);
00076 static int sinfo_utl_illumcorr(cpl_parameterlist *config, cpl_frameset *set);
00077 
00078 
00079 #define SINFO_DOUBLE_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
00080 
00081 static cpl_error_code 
00082 sinfo_tools_sort_double(
00083         double  *   pix_arr,
00084         int         n);
00085 
00086 static cpl_frame* 
00087 sinfo_get_dummy_object(cpl_frameset* obj_set);
00088 
00089 static void 
00090 sinfo_illumcorr_config_add (cpl_parameterlist *list);
00091 
00092 static int 
00093 create_illumcorr (const char* plugin_id, 
00094                  cpl_parameterlist *cpl_cfg, 
00095                  cpl_frameset* sof,
00096                  const char *name_i);
00097 static int 
00098 sinfo_illumcorr_create_bins (cpl_imagelist *sky, 
00099                  int llx, int lly, int urx, int ury,
00100                  int spec_bin, 
00101                  double min_flux,
00102                  int ** start,
00103                  int ** end,
00104                  int z1, int z2);
00105 
00106 static int 
00107 sinfo_juha_function1d_natural_spline(double *, double *, int, double *, 
00108                                      double *, int);
00109 
00110 static int 
00111 sinfo_function1d_search_value(double *, int, double, int *) ;
00112 
00113 static cpl_vector * 
00114 sinfo_vector_filter_median_create(const cpl_vector * v, int  hw);
00115 
00116 static cpl_vector * 
00117 sinfo_juha_vector_filter_median_create(const cpl_vector * v, int  hw);
00118 
00119 static double 
00120 sinfo_image_get_median_window  (const cpl_image *image, 
00121                 int llx, int lly, int urx, int ury);
00122 
00123 
00124 /*-----------------------------------------------------------------------------
00125                             Static variables
00126 -----------------------------------------------------------------------------*/
00127 
00128 static char sinfo_utl_illumcorr_description1[] =
00129 "This recipe calculates illumination correction based on sky emission.\n"
00130 "The input files are sky (or object) frames tagged\n"
00131 " SKY_NODDING (OBJECT_NODDING)\n"
00132 "Master calibration frames:\n";
00133 
00134 
00135 static char sinfo_utl_illumcorr_description2[] =
00136 "A corresponding (DIT) dark frame (tag=MASTER_DARK)"
00137 "A corresponding (band,preoptics) wavelength map image (tag=WAVE_MAP)\n"
00138 "A corresponding (band,preoptics) master flat field (tag=MASTER_FLAT_LAMP)\n"
00139 "A corresponding (band,preoptics) master bad pixel map (tag=MASTER_BP_MAP)\n"
00140 "A corresponding (band,preoptics) slitlets position frame (tag=SLIT_POS)\n"
00141 "A corresponding (band) distortion table (tag=DISTORTION)\n"
00142 "A corresponding (band) slitlet distance table (tag=SLITLETS_DISTANCE)\n";
00143 
00144 
00145 static char sinfo_utl_illumcorr_description3[] =
00146 "The output is a cube resulting from the analysis of sky emission\n";
00147 
00148 
00149 static char sinfo_utl_illumcorr_description4[] =
00150 "Information on relevant parameters can be found with\n"
00151 "esorex --params sinfo_utl_illumcorr\n"
00152 "esorex --help sinfo_utl_illumcorr\n"
00153 "\n";
00154 
00155 static char sinfo_utl_illumcorr_description[1300];
00156 
00157 /*-----------------------------------------------------------------------------
00158                                 Functions code
00159 -----------------------------------------------------------------------------*/
00160 
00161 /*---------------------------------------------------------------------------*/
00166 /*---------------------------------------------------------------------------*/
00168 static int
00169 sinfo_utl_illumcorr_create(cpl_plugin *plugin)
00170 {
00171 
00172   /*
00173    * We have to provide the option we accept to the application.
00174    * We need to setup our parameter list and hook it into the recipe
00175    * interface.
00176    */
00177   cpl_recipe *recipe = (cpl_recipe *)plugin;
00178   recipe->parameters = cpl_parameterlist_new();
00179   if(recipe->parameters == NULL) {
00180     return 1;
00181   }
00182   cpl_error_reset();
00183   irplib_reset();
00184 
00185   /*
00186    * Fill the parameter list.
00187    */
00188   sinfo_product_config_add (recipe->parameters);
00189   sinfo_prepare_stacked_frames_config_add(recipe->parameters);
00190   sinfo_objnod_config_add(recipe->parameters); 
00191   sinfo_illumcorr_config_add (recipe->parameters);
00192  
00193   return 0;
00194 
00195 }
00196 
00197 static int
00198 sinfo_utl_illumcorr_exec(cpl_plugin *plugin)
00199 {
00200  
00201   cpl_recipe *recipe = (cpl_recipe *) plugin;
00202   if(recipe->parameters == NULL) {
00203     return 1;
00204   }
00205   if(recipe->frames == NULL) {
00206     return 1;
00207   }
00208 
00209   return sinfo_utl_illumcorr(recipe->parameters, recipe->frames);
00210 
00211 }
00212 
00213 static int
00214 sinfo_utl_illumcorr_destroy(cpl_plugin *plugin)
00215 {
00216   cpl_recipe *recipe = (cpl_recipe *) plugin;
00217   /*
00218    * We just destroy what was created during the plugin initializzation phase
00219    * i.e. the parameter list. The frame set is managed by the application which
00220    * called us, so that we must not touch it.
00221    */
00222 
00223   cpl_parameterlist_delete(recipe->parameters);
00224 
00225   return 0;
00226 
00227 }
00228 
00229 
00230 int
00231 cpl_plugin_get_info(cpl_pluginlist *list)
00232 {
00233 
00234   cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe);
00235   cpl_plugin *plugin = &recipe->interface;
00236 
00237   strcpy(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description1);
00238   strcat(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description2);
00239   strcat(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description3);
00240   strcat(sinfo_utl_illumcorr_description,sinfo_utl_illumcorr_description4);
00241 
00242   cpl_plugin_init(plugin,
00243           CPL_PLUGIN_API,
00244           SINFONI_BINARY_VERSION,
00245           CPL_PLUGIN_TYPE_RECIPE,
00246           "sinfo_utl_illumcorr",
00247           "Object data reduction",
00248            sinfo_utl_illumcorr_description,
00249           "Juha Reunanen",
00250           "reunanen@strw.leidenuniv.nl",
00251           sinfo_get_license(),
00252           sinfo_utl_illumcorr_create,
00253           sinfo_utl_illumcorr_exec,
00254           sinfo_utl_illumcorr_destroy);
00255 
00256   cpl_pluginlist_append(list, plugin);
00257 
00258   return 0;
00259 
00260 }
00261 
00262 /*
00263  * The actual recipe actually start here.
00264  */
00265 
00266 static int
00267 sinfo_utl_illumcorr(cpl_parameterlist *config, cpl_frameset *set)
00268 {
00269   char outname[FILE_NAME_SZ];
00270 
00271   int i=0;
00272   int k=0;
00273 
00274   int ind=0;
00275   int nsky=0;
00276   int nobj=0;
00277   int ncdb=0;
00278   int nstk=0;
00279 
00280   cpl_frameset  * obj_set=NULL;
00281   cpl_frameset  * sky_set=NULL;
00282   cpl_frameset  * cdb_set=NULL;
00283   cpl_frameset  * wrk_set=NULL;
00284   cpl_frameset  * stk_set=NULL;
00285   cpl_frame     * sky_frm=NULL;
00286 
00287   cpl_frame     * dup_frm=NULL; 
00288   cpl_frame     * cdb_frm=NULL;
00289   cpl_frame     * wrk_frm=NULL;
00290   cpl_frameset  * ref_set=NULL;
00291 
00292   cpl_frame     * dark_frm=NULL;
00293 
00294   fake* fk;
00295 
00296 
00297   cpl_image           *   ima1=NULL ;
00298   cpl_image           *   ima2=NULL ;
00299   cpl_image           *   resima=NULL ;
00300   cpl_propertylist    *   plist=NULL ;
00301   cpl_frame           *   product_frame=NULL;
00302   const char *name_i=NULL;
00303 
00304   /* cpl_parameterlist_dump(config); */
00305         sinfo_msg("Welcome to SINFONI Pipeline release %d.%d.%d",
00306         SINFONI_MAJOR_VERSION,SINFONI_MINOR_VERSION,SINFONI_MICRO_VERSION);
00307 
00308   if(sinfo_dfs_set_groups(set)) {
00309     sinfo_msg_error("Cannot identify RAW and CALIB frames") ;
00310     return -1;
00311   }
00312 
00313   dark_frm = cpl_frameset_find(set,PRO_MASTER_DARK);
00314   if (dark_frm == NULL) {
00315     sinfo_msg_error("Cannot find dark frame") ;
00316     return (-1);
00317   }
00318 
00319   ref_set=cpl_frameset_duplicate(set);
00320 
00321   obj_set=cpl_frameset_new();
00322   sky_set=cpl_frameset_new();
00323   cdb_set=cpl_frameset_new();
00324   fk = sinfo_fake_new();
00325 
00326   sinfo_extract_obj_frames(set,obj_set);
00327   sinfo_extract_sky_frames(set,sky_set);
00328   sinfo_extract_mst_frames(set,cdb_set);
00329 
00330   nobj=cpl_frameset_get_size(obj_set);
00331   nsky=cpl_frameset_get_size(sky_set);
00332   ncdb=cpl_frameset_get_size(cdb_set);
00333 
00334   if ((nobj==0) && (nsky==0)) {
00335     sinfo_msg_error("Empty input set");
00336     cpl_frameset_delete(obj_set);
00337     cpl_frameset_delete(sky_set);
00338     cpl_frameset_delete(cdb_set);
00339     cpl_frameset_delete(ref_set);
00340     sinfo_fake_delete(&fk);
00341     return (-1);
00342   }
00343 
00344 
00345   /*
00346    *  Create median collapsed sky frame either from real SKY frames, 
00347    *  or from jittered OBJECT frames
00348    */
00349   if ( nsky != 0) {
00350     if( (sky_frm = sinfo_get_dummy_object(sky_set)) == NULL) {
00351       sinfo_msg_error("Problem to get dummy frame");
00352       cpl_frameset_delete(obj_set);
00353       cpl_frameset_delete(sky_set);
00354       cpl_frameset_delete(cdb_set);
00355       cpl_frameset_delete(ref_set);
00356       sinfo_fake_delete(&fk);
00357       return (-1);
00358     }
00359   }
00360   else {
00361     if( (sky_frm = sinfo_get_dummy_object(obj_set)) == NULL) {
00362       sinfo_msg_error("Problem to get dummy frame");
00363       cpl_frameset_delete(obj_set);
00364       cpl_frameset_delete(sky_set);
00365       cpl_frameset_delete(cdb_set);
00366       cpl_frameset_delete(ref_set);
00367       sinfo_fake_delete(&fk);
00368       return (-1);
00369     }
00370   }
00371   
00372   /* 
00373    *  Seems it's not possible to use draks as sky (due to INS.GRAT1.ENC)
00374    *  and stacking phase subtracts dark only in special circumstances...
00375    */
00376   ima1 = cpl_image_load(cpl_frame_get_filename(sky_frm),CPL_TYPE_FLOAT,0,0);
00377   ima2 = cpl_image_load(cpl_frame_get_filename(dark_frm),CPL_TYPE_FLOAT,0,0);
00378   resima = cpl_image_subtract_create(ima1, ima2);
00379   plist=cpl_propertylist_load(cpl_frame_get_filename(sky_frm), 0);
00380   cpl_image_delete(ima1);
00381   cpl_image_delete(ima2);
00382 
00383   product_frame = cpl_frame_new();
00384   cpl_frame_set_filename(product_frame, "out_fake_object2.fits") ;
00385   cpl_frame_set_tag(product_frame, "OBJECT_NODDING") ;
00386   cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
00387   cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_RAW) ;
00388   //cpl_frame_set_level(product_frame, CPL_FR) ;
00389   cpl_propertylist_erase_regexp(plist, "^ESO PRO CATG",0);
00390 
00391   cpl_image_save(resima, "out_fake_object2.fits", CPL_BPP_DEFAULT, plist,
00392          CPL_IO_DEFAULT) ;
00393   cpl_propertylist_delete(plist) ;
00394   cpl_image_delete(resima) ;
00395 
00396 
00397   /*
00398    *  Stack it - with some trickery...
00399    */ 
00400   wrk_set=cpl_frameset_new();
00401 
00402   dup_frm=cpl_frame_duplicate(product_frame);
00403   cpl_frame_set_tag   (dup_frm, "OBJECT_NODDING");
00404   cpl_frame_set_group (dup_frm ,CPL_FRAME_GROUP_RAW);
00405   cpl_frameset_insert(wrk_set,dup_frm);
00406 
00407   /* merge CDB frames to work set */
00408   for(k=0;k<ncdb;k++) {
00409     cdb_frm=cpl_frameset_get_frame(cdb_set,k);
00410     dup_frm=cpl_frame_duplicate(cdb_frm);
00411     cpl_frameset_insert(wrk_set,dup_frm);
00412   }
00413   
00414   
00415   /* defines a new name for the output stacked frame */ 
00416   sprintf(outname,"%s%d%s","out_stack",i,".fits");
00417   if(-1 == sinfo_new_stack_frames(config,wrk_set,
00418             PRO_OBJECT_NODDING_STACKED,i,fk,cpl_func)) {
00419     
00420     cpl_frameset_delete(wrk_set);
00421     //cpl_frameset_delete(tot_set);
00422     cpl_frameset_delete(obj_set);
00423     cpl_frameset_delete(sky_set);
00424     cpl_frameset_delete(cdb_set);
00425     cpl_frameset_delete(ref_set);
00426     sinfo_fake_delete(&fk);
00427     return -1;
00428   }
00429 
00430   stk_set=cpl_frameset_new();
00431   sinfo_contains_frames_kind(wrk_set,stk_set,PRO_STACKED);
00432   nstk=cpl_frameset_get_size(stk_set);
00433 
00434   for(k=0;k<nstk;k++) {
00435     wrk_frm=cpl_frameset_get_frame(stk_set,k);
00436     dup_frm = cpl_frame_duplicate(wrk_frm);
00437     cpl_frameset_insert(set,dup_frm);
00438   }
00439   cpl_frameset_delete(stk_set);
00440   cpl_frameset_delete(wrk_set);
00441 
00442   sinfo_msg("------------------------------") ; 
00443   sinfo_msg("CREATING SKY CUBE");
00444   sinfo_msg("------------------------------") ; 
00445   
00446 
00447     if ( -1 == (ind = sinfo_new_objnod(cpl_func,config, set, PRO_COADD_OBJ ) ) ) {
00448     sinfo_msg_error("NODDING SCIENCE FRAMES no. %d\n", ind) ;
00449     cpl_frameset_delete(obj_set);
00450     cpl_frameset_delete(sky_set);
00451     cpl_frameset_delete(cdb_set);
00452     cpl_frameset_delete(ref_set);
00453     sinfo_fake_delete(&fk);
00454     
00455     return (-1);
00456   }
00457   sinfo_msg("------------------------------") ; 
00458   sinfo_msg("CREATED SKY CUBE");
00459   sinfo_msg("------------------------------") ; 
00460   
00461 
00462   stk_set=cpl_frameset_new();
00463   sinfo_contains_frames_kind(set, stk_set, PRO_OBS_OBJ);
00464   nstk=cpl_frameset_get_size(stk_set);
00465 
00466   wrk_frm=cpl_frameset_get_frame(stk_set,0);  
00467   name_i = cpl_frame_get_filename(wrk_frm);
00468   
00469 
00470   //cpl_frameset_delete(wrk_set);
00471   cpl_frameset_delete(obj_set);
00472   cpl_frameset_delete(sky_set);
00473   cpl_frameset_delete(cdb_set);
00474   cpl_frameset_delete(ref_set);
00475   sinfo_fake_delete(&fk);
00476   cpl_frame_delete(sky_frm);
00477   create_illumcorr (cpl_func, config, set, name_i);
00478 
00479   return (0);
00480 
00481 }
00482 
00483 
00484 static cpl_frame*
00485 sinfo_get_dummy_object(cpl_frameset* obj_set)
00486 {
00487 
00488   cpl_imagelist* obj_list=NULL;
00489   cpl_image* fake_object=NULL;
00490   char filename[FILE_NAME_SZ];
00491   cpl_frame* frame=NULL;
00492   cpl_frame* object_frame=NULL;
00493 
00494   cpl_propertylist* plist=NULL;
00495  
00496   obj_list = cpl_imagelist_load_frameset(obj_set,CPL_TYPE_FLOAT,0,0);
00497   fake_object = cpl_imagelist_collapse_median_create(obj_list);
00498 
00499   frame = cpl_frameset_get_frame(obj_set,0);
00500   strcpy(filename,cpl_frame_get_filename(frame));
00501  
00502   if ((cpl_error_code)((plist = cpl_propertylist_load(filename, 0)) == NULL)) {
00503     sinfo_msg_error("getting header from reference ima frame %s",filename);
00504     cpl_propertylist_delete(plist) ;
00505     return NULL ;
00506   }
00507 
00508   if (cpl_propertylist_contains(plist, KEY_NAME_DPR_TYPE)) {
00509     cpl_propertylist_set_string(plist, KEY_NAME_DPR_TYPE, "OBJECT");
00510   } else {
00511     cpl_propertylist_append_string(plist, KEY_NAME_DPR_TYPE,"OBJECT") ;
00512   }
00513        
00514   if (cpl_image_save(fake_object, "out_fake_object.fits", CPL_BPP_DEFAULT, 
00515                    plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00516     sinfo_msg_error("Cannot save the product %s","out_fake_object.fits");
00517     cpl_propertylist_delete(plist) ;
00518     return NULL ;
00519   }
00520   cpl_propertylist_delete(plist);
00521        
00522   object_frame = cpl_frame_new() ;
00523   cpl_frame_set_filename(object_frame, "out_fake_object.fits") ;
00524   cpl_frame_set_tag(object_frame, "OBJECT") ;
00525   cpl_frame_set_type(object_frame, CPL_FRAME_TYPE_IMAGE);
00526   /*
00527   cpl_frame_set_group(object_frame, CPL_FRAME_GROUP_PRODUCT);
00528   */
00529   cpl_frame_set_level(object_frame, CPL_FRAME_LEVEL_FINAL);
00530   cpl_image_delete(fake_object);
00531   cpl_imagelist_delete(obj_list);
00532 
00533   return object_frame;
00534 }
00535 
00536 static void 
00537 sinfo_illumcorr_config_add (cpl_parameterlist *list) 
00538 {
00539   cpl_parameter *p;
00540 
00541   if (!list) {
00542     return;
00543   }
00544 
00545   p = cpl_parameter_new_range("sinfoni.illumcorr.spec_bin",
00546                               CPL_TYPE_INT,
00547                               "Number of spectral planes to be combined ",
00548                               "sinfoni.illumcorr",
00549                               100, 1, 200);
00550   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-spec_bin");
00551   cpl_parameterlist_append(list, p);
00552 
00553   p = cpl_parameter_new_value("sinfoni.illumcorr.min_flux",
00554                               CPL_TYPE_DOUBLE,
00555                               "Minimum flux in each spectral bin ",
00556                               "sinfoni.illumcorr",
00557                               0.0);
00558   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-min_flux");
00559   cpl_parameterlist_append(list, p);
00560 
00561   p = cpl_parameter_new_value("sinfoni.illumcorr.center_bins",
00562                               CPL_TYPE_BOOL,
00563                               "Center the spectral bins at prominent "
00564                               "emission features ",
00565                               "sinfoni.illumcorr",
00566                               FALSE);
00567   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-center_bins");
00568   cpl_parameterlist_append(list, p);
00569 
00570   p = cpl_parameter_new_enum("sinfoni.illumcorr.order",
00571                  CPL_TYPE_INT,
00572                  "The order of the polynomial to be fitted "
00573                              "for each slitlet",
00574                  "sinfoni.illumcorr",
00575                  0,
00576                  2,0,1);
00577   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-order");
00578   cpl_parameterlist_append(list, p);
00579   
00580   p = cpl_parameter_new_value("sinfoni.illumcorr.sigma",
00581                               CPL_TYPE_DOUBLE,
00582                               "Reject n-sigma deviant pixels on each slitlet ",
00583                               "sinfoni.illumcorr",
00584                               3.0);
00585   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-sigma");
00586   cpl_parameterlist_append(list, p);
00587 
00588   p = cpl_parameter_new_value("sinfoni.illumcorr.iterations",
00589                               CPL_TYPE_INT,
00590                               "Number of sigma rejection iterations to run ",
00591                               "sinfoni.illumcorr",
00592                               3);
00593   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-iter");
00594   cpl_parameterlist_append(list, p);
00595 
00596   p = cpl_parameter_new_range("sinfoni.illumcorr.llx",
00597                               CPL_TYPE_INT,
00598                               "Reference region coordinates ",
00599                               "sinfoni.illumcorr",
00600                               8, 0, 63);
00601   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-llx");
00602   cpl_parameterlist_append(list, p);
00603 
00604   p = cpl_parameter_new_range("sinfoni.illumcorr.lly",
00605                               CPL_TYPE_INT,
00606                               "Reference region coordinates ",
00607                               "sinfoni.illumcorr",
00608                               33, 0, 63);
00609   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-lly");
00610   cpl_parameterlist_append(list, p);
00611 
00612   p = cpl_parameter_new_range("sinfoni.illumcorr.urx",
00613                               CPL_TYPE_INT,
00614                               "Reference region coordinates ",
00615                               "sinfoni.illumcorr",
00616                               60, 0, 63);
00617   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-urx");
00618   cpl_parameterlist_append(list, p);
00619 
00620   p = cpl_parameter_new_range("sinfoni.illumcorr.ury",
00621                               CPL_TYPE_INT,
00622                               "Reference region coordinates ",
00623                               "sinfoni.illumcorr",
00624                               36, 0, 63);
00625   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-ury");
00626   cpl_parameterlist_append(list, p);
00627 
00628   p = cpl_parameter_new_enum("sinfoni.illumcorr.smooth0",
00629                  CPL_TYPE_INT,
00630                  "Smooth zeroth order terms by fitting "
00631                              "with polynomial (1),"
00632                  "with median filter (2) or not (0) ",
00633                  "sinfoni.illumcorr",
00634                  0, 
00635                  3, 0, 1, 2);
00636   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth0");
00637   cpl_parameterlist_append(list, p);
00638 
00639   p = cpl_parameter_new_value("sinfoni.illumcorr.smooth0_order",
00640                               CPL_TYPE_INT,
00641                               "Order of the smoothing polynomial for 0th term",
00642                               "sinfoni.illumcorr",
00643                               2);
00644   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth_order0");
00645   cpl_parameterlist_append(list, p);
00646 
00647   p = cpl_parameter_new_range("sinfoni.illumcorr.smooth0_size",
00648                               CPL_TYPE_INT,
00649                               "Size of the median filter for 0th "
00650                               "order smoothing ",
00651                               "sinfoni.illumcorr",
00652                               51,3, 301);
00653   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth0_size");
00654   cpl_parameterlist_append(list, p);
00655 
00656   p = cpl_parameter_new_value("sinfoni.illumcorr.smooth1",
00657                               CPL_TYPE_BOOL,
00658                               "Smooth higher (>0) order polynomials ",
00659                               "sinfoni.illumcorr",
00660                               TRUE);
00661   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth");
00662   cpl_parameterlist_append(list, p);
00663 
00664   p = cpl_parameter_new_value("sinfoni.illumcorr.smooth1_order",
00665                               CPL_TYPE_INT,
00666                               "Smoothing order for higher terms ",
00667                               "sinfoni.illumcorr",
00668                               2);
00669   cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,"illumcorr-smooth_order");
00670   cpl_parameterlist_append(list, p);
00671 
00672   p = cpl_parameter_new_value("sinfoni.illumcorr.illumcorr_sigma",
00673                               CPL_TYPE_DOUBLE,
00674                               "Reject all fits for which the rms is "
00675                               "illumcorr-sigma times bigger than the "
00676                               "median rms in each spectral bin" ,
00677                               "sinfoni.illumcorr",
00678                               5.0);
00679   cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,
00680                             "illumcorr-illumcorr_sigma");
00681   cpl_parameterlist_append(list, p);
00682 
00683 }
00684 
00685 static int
00686 create_illumcorr (const char* plugin_id,
00687           cpl_parameterlist *cpl_cfg,
00688           cpl_frameset* sof,
00689           const char *name_i)
00690 {
00691   cpl_parameter *p=NULL;
00692   double        min_flux=0;
00693   double        sigma=0;
00694   double        il_sigma=0;
00695   int           spec_bin=0;
00696   int           _order=0;
00697   cpl_imagelist *sky=NULL;
00698   cpl_imagelist *binnedsky=NULL;
00699   cpl_imagelist *result=NULL;
00700   cpl_image     *temp_image=NULL;
00701   cpl_image     *temp_image2=NULL;
00702   int            nplanes=0;
00703   int            i=0;
00704   int            j=0;
00705   int            k=0;
00706   int            kk=0;
00707   int            n=0;
00708   int            slitlet=0;
00709   int            bin=0;
00710   double        *median=NULL;
00711   double        *pos=NULL;
00712   double         temp=0;
00713   double         temp2=0;
00714   double        *inter_pos=NULL;
00715   double        *inter_val=NULL;
00716   double        *plane_pos=NULL;
00717   double        *plane_val=NULL;
00718   int            llx=0;
00719   int            lly=0;
00720   int            urx=0;
00721   int            ury=0;
00722   int            smooth_order=0;
00723   int            iter=0;
00724   cpl_vector    *row=NULL;
00725   cpl_vector    *model=NULL;
00726   cpl_vector    *xpos=NULL;
00727   cpl_vector    *tempvector=NULL;
00728   cpl_vector *tempvector2=NULL;
00729   double         mse=0.0;
00730   double         stddev=0.0;
00731   cpl_polynomial*poly=NULL;
00732   cpl_polynomial *poly2=NULL;
00733   double        *temparray=NULL;
00734   double  *tempxarray=NULL;
00735   double  * tempsarray=NULL;
00736   cpl_polynomial**coeffs=NULL;
00737   float         *data=NULL;
00738   double        *rms_values=NULL;
00739   double         rms_array[32];
00740   int smooth=0;
00741   int smooth0=0;
00742   int smooth_order0=0;
00743   int smooth_size0=0;
00744   int center_bins = 0;
00745 
00746   int *bin_start=NULL;
00747   int *bin_end=NULL;
00748   int z1=0;
00749   int z2=0;
00750   int nbins=0;
00751 
00752   FILE          *dumpfile=NULL;
00753 
00754   int order[32];
00755   int ok[64];
00756   int nbad=0;
00757 
00758 
00759   /*
00760    *  Get parameters
00761    */
00762   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.spec_bin");
00763   spec_bin = cpl_parameter_get_int(p);
00764   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.min_flux");
00765   min_flux =  cpl_parameter_get_double(p);
00766   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.order");
00767   _order = cpl_parameter_get_int(p);
00768   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.sigma");
00769   sigma = cpl_parameter_get_double(p);
00770   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.llx");
00771   llx = cpl_parameter_get_int(p);
00772   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.lly");
00773   lly = cpl_parameter_get_int(p);
00774   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.urx");
00775   urx = cpl_parameter_get_int(p);
00776   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.ury");
00777   ury = cpl_parameter_get_int(p);
00778   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.illumcorr_sigma");
00779   il_sigma = cpl_parameter_get_double(p);
00780 
00781   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth0");
00782   smooth0 = cpl_parameter_get_int (p);
00783   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth0_order");
00784   smooth_order0 = cpl_parameter_get_int (p);
00785   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth0_size");
00786   smooth_size0 = cpl_parameter_get_int (p);
00787 
00788   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth1");
00789   smooth = cpl_parameter_get_bool (p);
00790   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.smooth1_order");
00791   smooth_order = cpl_parameter_get_int (p);
00792 
00793   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.iterations");
00794   iter = cpl_parameter_get_int (p);
00795 
00796   p = cpl_parameterlist_find(cpl_cfg, "sinfoni.illumcorr.center_bins");
00797   center_bins = cpl_parameter_get_bool (p);
00798 
00799 /*   cpl_msg_set_level(CPL_MSG_DEBUG); */
00800 
00801   /*
00802    *  Allocate resources
00803    */
00804   sky       = cpl_imagelist_load(name_i, CPL_TYPE_FLOAT, 0);
00805   nplanes   = cpl_imagelist_get_size(sky);
00806 
00807   /*  Determine the start and end points of data within the 
00808    *  reference region  */   
00809   z1 = 0;
00810   z2=nplanes -1;
00811   while (z1<nplanes 
00812      && isnan(cpl_image_get_mean_window(cpl_imagelist_get(sky, z1),
00813                         llx, lly, urx, ury)))
00814     z1++;
00815   while (z2>=0
00816      && isnan(cpl_image_get_mean_window(cpl_imagelist_get(sky, z2),
00817                         llx, lly, urx, ury)))
00818     z2--;
00819   z1 += 2;
00820   z2 -= 2;
00821   if (z1>=nplanes || z2 <0 || z2<=z1) {
00822     sinfo_msg_error ("Start z = %d, end z = %d", z1, z2);
00823     cpl_imagelist_delete (sky);
00824     return (-1);
00825   }
00826   sinfo_msg ("Start z = %d, end z = %d", z1, z2);
00827 
00828   binnedsky = cpl_imagelist_new ();
00829   median    = (double*) cpl_calloc(nplanes/spec_bin, sizeof(double));
00830   pos       = (double*) cpl_calloc(nplanes/spec_bin, sizeof(double));
00831   temparray = (double*) cpl_calloc(64, sizeof(double));
00832   tempxarray= (double*) cpl_calloc(64, sizeof(double));
00833   tempsarray= (double*) cpl_calloc (nplanes/spec_bin, sizeof(double));
00834   plane_pos = (double*) cpl_calloc (nplanes/spec_bin, sizeof(double));
00835   plane_val = (double*) cpl_calloc (nplanes/spec_bin, sizeof(double));
00836   coeffs    = (cpl_polynomial**) cpl_calloc(32*(nplanes/spec_bin), 
00837                         sizeof(cpl_polynomial*));
00838   rms_values= (double*) cpl_calloc (32*(nplanes/spec_bin), sizeof (double));
00839   inter_pos = (double*) cpl_calloc (nplanes, sizeof(double));
00840   inter_val = (double*) cpl_calloc (nplanes, sizeof(double));
00841 
00842   model     = cpl_vector_new(64);
00843   xpos      = cpl_vector_new(64);
00844 
00845   for (i=0; i<64; i++)
00846     cpl_vector_set(xpos, i, (double)(i)-((double)urx-(double)llx)/2.0);
00847   for (i=0; i<nplanes; i++)
00848     inter_pos[i] = (double)i;
00849 
00850   /*
00851    *  This array could be given as input file for the recipe.
00852    *  Generally, 0th order fitting is sufficient (and of course
00853    *  more robust), but few slitlets might require 1st order.
00854    */
00855   for (i=0; i<32; i++)
00856     order[i] = _order;
00857 
00858 
00859   if (center_bins == 1) {
00860     sinfo_msg("Using centering on emission features\n");
00861     nbins = sinfo_illumcorr_create_bins (sky,llx, lly, urx, ury,
00862                      spec_bin, min_flux,
00863                      &bin_start, &bin_end,
00864                      z1, z2);
00865   }
00866   else {
00867     sinfo_msg("Using simple spectral binning - "
00868           "not centering on emission features\n");
00869     nbins = (z2-z1)/spec_bin;
00870     bin_start = (int*)cpl_calloc(nbins+1, sizeof(int));
00871     bin_end =  (int*)cpl_calloc(nbins+1, sizeof(int));
00872     for (i = 0; i<nbins; i++) {
00873       bin_start[i] = z1+i*spec_bin;
00874       bin_end[i]   = z1+(i+1)*spec_bin - 1;
00875     }
00876     if (bin_end[nbins-1]<z2-spec_bin/10) {
00877       bin_start[nbins] = bin_end[nbins-1]+1;
00878       bin_end[nbins] = z2;
00879       nbins++;
00880     }
00881   }
00882   
00883   /*
00884    *  - bin the cube in spectral direction
00885    *  - calculate the median (=reference value) in region 
00886    *    (llx,lly) - (urx,ury)
00887    *  - calculate the weighted position of the each spectral bin
00888    */
00889   sinfo_msg("Binning the cube and calculating statistics\n");
00890   for (i=0; i<nbins; i++) {
00891     temp_image = cpl_image_duplicate(cpl_imagelist_get(sky, bin_start[i]));
00892     median[i] = sinfo_image_get_median_window (temp_image, llx, lly, urx, ury);
00893     pos[i] = median[i] * (double)bin_start[i];
00894     cpl_imagelist_set (binnedsky, temp_image, i);
00895     for (j=bin_start[i]+1; j<bin_end[i]; j++) {
00896       temp_image2 = cpl_imagelist_get (sky, j);
00897       cpl_image_add (temp_image, temp_image2);
00898       temp = sinfo_image_get_median_window (temp_image2, llx, lly, urx, ury);
00899       median[i] = median[i] + temp;
00900       pos[i] = pos[i] + temp*(double)j;
00901     }
00902     temp2 =(double)(bin_end[i]-bin_start[i]+1);
00903     cpl_image_divide_scalar (temp_image, temp2);
00904     pos[i] = pos[i]/median[i];
00905     median[i] = median[i] / temp2;
00906     sinfo_msg_debug("median image=%g at %g",median[i], pos[i]); 
00907   }
00908   
00909   sinfo_msg("Fitting slitlets between x=%d - x=%d\n", llx, urx);
00910   sinfo_msg("Fitting order %d\n", _order); 
00911   for (k=0;k<nbins; k++) {
00912     if (median[k]>min_flux) {
00913       for (j=0; j<32; j++) {
00914     row=cpl_vector_new_from_image_row(cpl_imagelist_get(binnedsky,k),2*j+1);
00915     n = 0;
00916     for (i=llx; i<=urx; i++) {
00917       if (!isnan(cpl_vector_get(row, i))) {
00918         ok[i] = 1;
00919         temparray[n]  = cpl_vector_get(row, i);
00920         tempxarray[n] = cpl_vector_get(xpos, i);
00921         n++;
00922       }
00923       else
00924         ok[i] = 0;
00925     }
00926 
00927     /* The ends of cube are always filled with NaNs => n==0*/
00928     if (n>20) {
00929       tempvector = cpl_vector_wrap (n, temparray);
00930       tempvector2= cpl_vector_wrap (n, tempxarray);
00931       poly = cpl_polynomial_fit_1d_create (tempvector2, tempvector, 
00932                                                order[j], &mse);
00933       
00934       if (poly == NULL)
00935         sinfo_msg("Fitting failed (plane %d, row %d): %s",
00936                  k, j, (char* ) cpl_error_get_message());
00937       else {
00938         if (sigma>0 && iter>0) {
00939           for (kk = 0; kk<iter; kk++) {
00940         cpl_vector_fill_polynomial (model, poly, 0.0, 1.0);
00941         cpl_vector_subtract (model, row);
00942         
00943         /* Calculate stdev NaN-correctly */
00944         n = 0;
00945         for (i=llx; i<=urx; i++) 
00946           if (ok[i] == 1) 
00947             temparray[n++]  = cpl_vector_get(model, i);
00948         stddev = cpl_vector_get_stdev(tempvector);
00949         
00950         for (i=llx; i<=urx; i++) 
00951           if (ok[i] == 1) 
00952             if (fabs(cpl_vector_get(model, i))>(stddev*sigma))
00953               ok[i] = 0;
00954         
00955         
00956         n = 0;
00957         for (i=llx; i<=urx; i++) {
00958           if (ok[i] == 1) {
00959             temparray[n]  = cpl_vector_get(row, i);
00960             tempxarray[n] = cpl_vector_get(xpos, i);
00961             n++;
00962           }
00963         }
00964         cpl_polynomial_delete(poly);
00965         if (n>20) {
00966           cpl_vector_unwrap (tempvector);
00967           cpl_vector_unwrap (tempvector2);
00968           tempvector = cpl_vector_wrap (n, temparray);
00969           tempvector2= cpl_vector_wrap (n, tempxarray);
00970           stddev = cpl_vector_get_stdev(tempvector);
00971           
00972           poly = cpl_polynomial_fit_1d_create (tempvector2, 
00973                                                        tempvector, 
00974                                                        order[j], &mse);
00975           if (poly == NULL)
00976             break;
00977         }
00978         else {
00979           poly = NULL;
00980           break;
00981         }
00982 /*        printf ("%d %e ", n, stddev); */
00983           }
00984         }
00985 
00986         if (poly!=NULL) {
00987           coeffs[j*nbins+k] = poly;
00988           rms_values[j*nbins+k] = sqrt(mse);
00989         }
00990         else
00991           coeffs[j*nbins+k] = NULL;
00992       }
00993       cpl_vector_unwrap (tempvector);
00994       cpl_vector_unwrap (tempvector2);
00995     }
00996     cpl_vector_delete(row);
00997       }
00998     }
00999   }
01000 
01001   /*
01002    *   These should (probably) be saved in a fits file...
01003    */
01004   sinfo_msg("Writing dat out_illum.dat\n");
01005   dumpfile = fopen ("out_illum.dat","w");
01006   fprintf (dumpfile, "# slitlet, pos, median, rms, coeff0, coeff1...\n");
01007   for (slitlet = 0; slitlet<32; slitlet++)
01008     for (bin=0; bin<nbins; bin++) {
01009       poly = coeffs[slitlet*nbins+bin];
01010       if (poly != NULL) {
01011     fprintf (dumpfile, "%d %f %f %f ",slitlet, pos[bin],
01012          median[bin],
01013          rms_values[slitlet*nbins+bin]);
01014     for (i=0; i<(cpl_polynomial_get_degree(poly)+1); i++) {
01015       temp = cpl_polynomial_get_coeff (poly, &i);
01016       fprintf (dumpfile, "%f ", temp);
01017     }
01018     fprintf (dumpfile, "\n");
01019       }
01020     }
01021   fclose (dumpfile);
01022 
01023   /*
01024    *    Remove poor fits:
01025    *      - calculate the median rms of all fits
01026    *      - throw away the fits whose rms is il_sigma*median_rms
01027    */
01028   sinfo_msg("Removing poor fits - factor %f", il_sigma);
01029   n = 0;
01030   i = 0;
01031   nbad=0;
01032   sinfo_msg("max loop over bin =%d",nbins);
01033   for (bin=0; bin<nbins; bin++) {
01034     k = 0;
01035     for (slitlet=0; slitlet<32; slitlet++)
01036       if (coeffs[slitlet*nbins+bin] != NULL)
01037     rms_array[k++] = rms_values[slitlet*nbins+bin];
01038     if (k>0) {
01039       /* For some bizarre reason, cpl_tools_get_median_double returns 
01040        *  -1076245448.000000 (is that NaN?). On closer inspection,
01041        * it seems to have replaced one of the numbers in array with NaN...*/
01042       tempvector = cpl_vector_wrap (k, &rms_array[0]);
01043       temp = cpl_vector_get_median (tempvector);
01044       sinfo_msg("median temp=%g",temp);
01045       cpl_vector_unwrap (tempvector);
01046       for (slitlet=0; slitlet<32; slitlet++) 
01047     if (coeffs[slitlet*nbins+bin] != NULL) {
01048       i++;
01049       if (rms_values[slitlet*nbins+bin]>(il_sigma*temp)) {
01050         cpl_polynomial_delete(coeffs[slitlet*nbins+bin]);
01051         coeffs[slitlet*nbins+bin] = NULL;
01052         n++;
01053       }
01054     } else {
01055           nbad++;
01056     }
01057         
01058     }
01059   }
01060   sinfo_msg("Removed %d poor fits out of %d. Bad coeffs=%d", n,i,nbad);
01061 
01062   if(smooth0 == 1) {
01063     sinfo_msg("Smoothing zeroth terms (order %d)", smooth_order0);
01064     /*
01065      *  Since the new centering scheme will pro
01066      */ 
01067     for (slitlet = 0; slitlet<32; slitlet++) {
01068       k = 0;
01069       for (bin=0; bin<nbins; bin++) {
01070     if (coeffs[slitlet*nbins+bin] != NULL) {
01071       poly = coeffs[slitlet*nbins+bin];
01072       i = 0;
01073       temp = cpl_polynomial_get_coeff (poly, &i);
01074       plane_pos[k] = pos[bin];
01075       plane_val[k] = temp/median[bin];
01076       k++;
01077     }
01078       }
01079       if (k>0) {
01080     tempvector = cpl_vector_wrap (k, plane_pos);
01081     tempvector2= cpl_vector_wrap (k, plane_val);
01082     poly2 = cpl_polynomial_fit_1d_create (tempvector, tempvector2, 
01083                                               smooth_order0, &mse);
01084     cpl_vector_unwrap (tempvector);
01085     cpl_vector_unwrap (tempvector2);
01086     for (bin=0; bin<nbins; bin++) {
01087       if (coeffs[slitlet*nbins+bin] != NULL) {
01088         poly = coeffs[slitlet*nbins+bin];
01089         i = 0;
01090         temp2 = cpl_polynomial_eval_1d (poly2, pos[bin], NULL);
01091         cpl_polynomial_set_coeff (poly, &i, temp2*median[bin]);
01092       }
01093     }
01094     cpl_polynomial_delete(poly2); 
01095       }
01096       else 
01097     sinfo_msg_warning ("Not enough data points in slitlet %d", slitlet);
01098     }
01099   }
01100   else if (smooth0 == 2) {
01101     sinfo_msg("Smoothing zeroth terms (median filter size %d)", smooth_size0);
01102     smooth_size0 = smooth_size0/2;
01103     if (smooth_size0 <= 0) smooth_size0 = 1;
01104     for (slitlet = 0; slitlet<32; slitlet++) {
01105       k = 0;
01106       for (bin=0; bin<nbins; bin++) {
01107     if (coeffs[slitlet*nbins+bin] != NULL) {
01108       poly = coeffs[slitlet*nbins+bin];
01109       i = 0;
01110       temp = cpl_polynomial_get_coeff (poly, &i);
01111       //plane_pos[k] = pos[bin];
01112       plane_val[k] = temp/median[bin];
01113       k++;
01114     }
01115       }
01116       if (k>0) {
01117     tempvector = cpl_vector_wrap (k, plane_val);
01118     tempvector2= sinfo_juha_vector_filter_median_create (tempvector, 
01119                                                              smooth_size0);
01120     cpl_vector_unwrap (tempvector);
01121     kk = 0;
01122     for (bin=0; bin<nbins; bin++) {
01123       if (coeffs[slitlet*nbins+bin] != NULL) {
01124         poly = coeffs[slitlet*nbins+bin];
01125         i = 0;
01126         cpl_polynomial_set_coeff(poly, &i, cpl_vector_get(tempvector2, kk++)
01127                       *median[bin]);
01128       }
01129     }
01130     cpl_vector_delete (tempvector2);
01131       }
01132     }
01133   }
01134 
01135   if(smooth == 1) {
01136     sinfo_msg("Smoothing higher terms (with order %d)", smooth_order);
01137     for (slitlet = 0; slitlet<32; slitlet++) {
01138       if (order[slitlet]>0) {
01139     for (j=1; j<=order[slitlet]; j++) {
01140       k = 0;
01141       for (bin=0; bin<nbins; bin++) {
01142         if (coeffs[slitlet*nbins+bin] != NULL) {
01143           poly = coeffs[slitlet*nbins+bin];
01144           i = 0;
01145 /*        temp = cpl_polynomial_get_coeff (poly, &i); */
01146           temp2 = cpl_polynomial_get_coeff (poly, &j);
01147           plane_pos[k] = pos[bin];
01148           plane_val[k] = temp2/median[bin];
01149           k++;
01150         }
01151       }
01152       if (k>0) {
01153         tempvector = cpl_vector_wrap (k, plane_pos);
01154         tempvector2= cpl_vector_wrap (k, plane_val);
01155         poly2 = cpl_polynomial_fit_1d_create (tempvector, tempvector2, 
01156                                                   smooth_order, &mse);
01157         cpl_vector_unwrap (tempvector);
01158         cpl_vector_unwrap (tempvector2);
01159         for (bin=0; bin<nbins; bin++) {
01160           if (coeffs[slitlet*nbins+bin] != NULL) {
01161         poly = coeffs[slitlet*nbins+bin];
01162           i = 0;
01163           /*          temp = cpl_polynomial_get_coeff (poly, &i); */
01164           temp2 = cpl_polynomial_eval_1d (poly2, pos[bin], NULL);
01165           cpl_polynomial_set_coeff (poly, &j, temp2*median[bin]);
01166           }
01167         }
01168         cpl_polynomial_delete(poly2); 
01169       }
01170       else 
01171         sinfo_msg_warning ("Not enough data points in slitlet %d\n", 
01172                                slitlet);
01173     }
01174       }
01175     }
01176   }
01177 
01178   sinfo_msg("Creating cube for illumination correction\n");
01179   result = cpl_imagelist_new ();
01180   for (i=0; i<nplanes; i++) {
01181     temp_image = cpl_image_new (64, 64, CPL_TYPE_FLOAT);
01182     cpl_imagelist_set (result, temp_image, i);
01183   }
01184     
01185   sinfo_msg("ok1 nplanes=%d spec_bin=%d",nplanes,spec_bin);
01186   if ( nbins>5) {
01187     sinfo_msg("Interpolating\n");
01188     for (slitlet = 0; slitlet<32; slitlet++) {
01189       for (i=0; i<64; i++) {
01190     k = 0;
01191     for (bin=0; bin<nbins; bin++) {
01192       if (coeffs[slitlet*nbins+bin] != NULL) {
01193         plane_pos[k] = pos[bin];
01194         plane_val[k] = cpl_polynomial_eval_1d(coeffs[slitlet*nbins+bin], 
01195                           cpl_vector_get(xpos, i),NULL)/
01196                                                   median[bin];
01197         k++;
01198       }
01199     }
01200         
01201     if (k>3) {
01202       sinfo_juha_function1d_natural_spline (plane_pos, plane_val, k,
01203                          &inter_pos[(int)plane_pos[0]],
01204                          &inter_val[(int)plane_pos[0]],
01205                          (int)(plane_pos[k-1]-plane_pos[0]));
01206       for (j=0; j<=(int)plane_pos[0]; j++)
01207         inter_val[j] = inter_val[(int)plane_pos[0]+1];
01208       for (j=(int)plane_pos[k-1]-1; j<nplanes; j++)
01209         inter_val[j] = inter_val[(int)plane_pos[k-1]-2];
01210       for (k=0; k<nplanes; k++) {
01211         data = cpl_image_get_data_float(cpl_imagelist_get(result, k));
01212         data[i + (2*slitlet)*64] = inter_val[k];
01213         data[i + (2*slitlet+1)*64] = inter_val[k];
01214             /*sinfo_msg("inter_val=%g",inter_val[k]);*/
01215       }
01216     }
01217     else
01218     sinfo_msg_warning ("Too few points %d\n", slitlet);
01219       }
01220     }
01221   }
01222   else if (nbins==1) {
01223     sinfo_msg("Filling the illumination cube\n");
01224     for (slitlet = 0; slitlet<32; slitlet++) {
01225       for (i=0; i<64; i++) {
01226     if (coeffs[slitlet] != NULL) {
01227       temp = cpl_polynomial_eval_1d(coeffs[slitlet], 
01228                     cpl_vector_get(xpos, i),NULL)/median[0];
01229       for (k=0; k<nplanes; k++) {
01230         data = cpl_image_get_data_float(cpl_imagelist_get(result, k));
01231         data[i + (2*slitlet)*64] = temp;
01232         data[i + (2*slitlet+1)*64] = temp;
01233             sinfo_msg("temp=%g",temp);
01234       }
01235     }
01236       }
01237     }
01238   } else {
01239 
01240   }
01241   sinfo_msg("ok2");
01242 
01243 
01244   sinfo_msg("Writing ima out_illum.fits\n");
01245 /*   pl = cpl_propertylist_load (name_i, 0); */
01246 /*   if (cpl_propertylist_contains(pl, KEY_NAME_PRO_CATG)) */
01247 /*     cpl_propertylist_set_string (pl, KEY_NAME_PRO_CATG, PRO_ILL_COR); */
01248 /*   else */
01249 /*     cpl_propertylist_append_string (pl, KEY_NAME_PRO_CATG, PRO_ILL_COR); */
01250 
01251 /*   cpl_imagelist_save(result, "out_illum.fits", CPL_BPP_DEFAULT, pl, 0); */
01252 
01253   sinfo_pro_save_ims(result,sof,sof,"out_illum.fits", 
01254              PRO_ILL_COR,NULL,plugin_id, cpl_cfg); 
01255 
01256   /*
01257    *   These should (probably) be saved in a fits file...
01258    */
01259   sinfo_msg("Writing dat out_illum2.dat\n");
01260   dumpfile = fopen ("out_illum2.dat","w");
01261   fprintf (dumpfile, "# slitlet, pos, median, rms, coeff0, coeff1...\n");
01262   for (slitlet = 0; slitlet<32; slitlet++)
01263     for (bin=0; bin<nbins; bin++) {
01264       poly = coeffs[slitlet*nbins+bin];
01265       if (poly != NULL) {
01266     fprintf (dumpfile, "%d %f %f %f ",slitlet, pos[bin],
01267          median[bin],
01268          rms_values[slitlet*nbins+bin]);
01269     for (i=0; i<(cpl_polynomial_get_degree(poly)+1); i++) {
01270       temp = cpl_polynomial_get_coeff (poly, &i);
01271       fprintf (dumpfile, "%f ", temp);
01272     }
01273     fprintf (dumpfile, "\n");
01274       }
01275     }
01276   fclose (dumpfile);
01277 
01278   /*
01279    *  Clean up...
01280    */
01281   for (i = 0; i<32*nbins; i++)
01282     if (coeffs[i] != NULL)
01283       cpl_polynomial_delete(coeffs[i]);
01284   cpl_imagelist_delete (sky);
01285   cpl_imagelist_delete (binnedsky);
01286   cpl_imagelist_delete (result);
01287   cpl_free (pos);
01288   cpl_free (median);
01289   cpl_free (temparray);
01290   cpl_free (tempxarray);
01291   cpl_free (tempsarray);
01292   cpl_free (coeffs);
01293   cpl_free (inter_pos);
01294   cpl_free (inter_val);
01295   cpl_free (plane_pos);
01296   cpl_free (plane_val);
01297   cpl_free (rms_values);
01298   cpl_vector_delete (xpos);
01299   cpl_vector_delete (model);
01300 
01301   cpl_free (bin_start);
01302   cpl_free (bin_end);
01303 
01304   return (1);
01305 }
01306 
01307 /*
01308  *  sinfo_illumcorr_create_bins:
01309  *    - searches for the sky emission lines
01310  *    - increases the size of the bin to include two or more emission
01311  *      lines if they are too close to each other
01312  *    - fills the space between emission lines with bins if
01313  *      thermal background has enough flux
01314  *    - copies the start and end points of bins to two arrays
01315  *      (returned in **start and **end)
01316  *
01317  *  Returns: the number bins created
01318  *
01319  *  The arrays start and end must be deallocated with cpl_free()
01320  */
01321 static int
01322 sinfo_illumcorr_create_bins (cpl_imagelist *sky, 
01323                  int llx, int lly, int urx, int ury,
01324                  int spec_bin, 
01325                  double min_flux,
01326                  int ** start,
01327                  int ** end,
01328                  int z1, int z2)
01329 {
01330   int temp_i=0;
01331   double testarray3[15];
01332   double temp_double=0;
01333   int  i=0, j=0, k=0,kk=0,nplanes=0; 
01334 
01335   int norig = 0, nmerged = 0, ncont = 0, nline=0;
01336 
01337   int *pos=NULL;
01338   int *x1=NULL;
01339   int *x2=NULL;
01340   int *x1b=NULL;
01341   int *x2b=NULL;
01342   int *s1=NULL;
01343   int *s2=NULL;
01344   double *flux=NULL;
01345   double *spec=NULL;
01346   double  *spec_cont=NULL;
01347   double *spec_line=NULL;
01348 
01349 
01350   cpl_image     *temp_image=NULL;
01351 
01352   nplanes   = cpl_imagelist_get_size(sky);
01353 
01354   spec = (double*) cpl_calloc(nplanes, sizeof(double));
01355   spec_line = (double*) cpl_calloc(nplanes, sizeof(double));
01356   spec_cont = (double*) cpl_calloc(nplanes, sizeof(double));
01357 
01358   /*  there should be no way of actually needing this large arrays*/
01359   pos = (int*) cpl_calloc(nplanes, sizeof(int));
01360   flux = (double*) cpl_calloc(nplanes, sizeof(double));
01361   x1 = (int*) cpl_calloc(nplanes, sizeof(int));
01362   x2 = (int*) cpl_calloc(nplanes, sizeof(int));
01363   x1b = (int*) cpl_calloc(nplanes, sizeof(int));
01364   x2b = (int*) cpl_calloc(nplanes, sizeof(int));
01365 
01366   for (i=z1; i<=z2; i++) {
01367     temp_image =  cpl_imagelist_get(sky, i);
01368     spec[i] = sinfo_image_get_median_window (temp_image, llx, lly, urx, ury);
01369   }
01370   for (i=z1+7; i<=z2-7; i++) {
01371     k = 0;
01372     for (j=-7; j<=7; j++)
01373       if (!isnan(spec[i+j]))
01374     testarray3[k++] = spec[i+j];
01375     if (k>0) {
01376       sinfo_tools_sort_double (&testarray3[0], k);
01377       spec_cont[i] =  testarray3[1];
01378     }
01379     else
01380       spec_cont[i] =  0./0.;
01381   }
01382 
01383   sinfo_msg_debug("Calculating pure line flux at pos: "
01384                   "original, continuum, line");
01385   for (i=z1; i<=z2; i++) {
01386     spec_line[i] = spec[i] - spec_cont[i];
01387     sinfo_msg_debug("Flux at %i = %g %g %g",
01388                      i,spec[i],spec_cont[i], spec_line[i]);   
01389   }
01390 
01391 
01392   /*
01393    *  Search for peaks
01394    */  
01395   sinfo_msg ("Searching for peaks");
01396   temp_double = -10000.0;
01397   i = z1+2; 
01398   while (i<=z2-2) {
01399     if (!isnan (spec_line[i])) {
01400       if (temp_double<spec_line[i]) {
01401         temp_i = i;
01402         temp_double = spec_line[i];
01403       }
01404       else {
01405         /*  Found a peak! */ 
01406         if (temp_i == i-1 && spec_line[temp_i]>min_flux) {
01407       k = 0;
01408       for (j=-spec_bin/2; j<=spec_bin/2; j++) 
01409         if (j+i>=0 && i+j<nplanes && isnan(spec[i+j])) {
01410           k = 1;
01411           break;
01412         }
01413       if (k==0) {
01414         pos[norig] = temp_i; // - spec_bin/2;
01415         flux[norig]   = temp_double;
01416         x1[norig] = temp_i;
01417         x2[norig] = temp_i;
01418         temp_double = -10000.0; 
01419         while (spec_line[i]<spec_line[i-1])
01420           i++;
01421         i--;
01422         norig++;
01423       }
01424         }
01425       }
01426     }      
01427     i++;
01428   }
01429 
01430   /*
01431    *   Merge the features which are too close to each other
01432    */
01433   sinfo_msg ("Merging emission features too close to each other");
01434   i = 0;
01435   while (i<norig) {
01436     if (flux[i] > 0.0) {
01437       j = i+1;
01438       while (j<norig 
01439          && (x1[j]-x2[i]) <=spec_bin 
01440          && flux[j]>0.0) {
01441     flux[j] = -100.0;
01442     x2[i] = x1[j];
01443     j++;
01444     nmerged++;
01445       }
01446     }
01447     i++;
01448   }
01449 
01450   nline = norig - nmerged;
01451 
01452   j = 0;
01453   for (i=0; i<norig; i++)
01454     if (flux[i]>0.0) {
01455       x1b[j] = x1[i] - spec_bin/2;
01456       x2b[j] = x2[i] + spec_bin/2;
01457       j++;
01458       /* sinfo_msg ("Bin start: %i end %i", x1[i], x2[i]); */
01459     }
01460 
01461   x1b[j] = nplanes+1;
01462 
01463   /*
01464    *  Check whether there is enough continuum (thermal background)
01465    *  for binning
01466    */
01467   j=0;
01468   i=z1;
01469   while (i<=z2) {
01470     if (!isnan (spec[i])) {
01471       if (x1b[j]-i < spec_bin) {
01472     i = x2b[j]+1;
01473     j++;
01474       }
01475       else {
01476     kk = 0;
01477     for (k=0; k<spec_bin; k++) 
01478       if (spec[i+k]>min_flux)
01479         kk++;
01480     if (kk==spec_bin) {
01481       x1[ncont] = i;
01482       x2[ncont] = i+spec_bin-1;
01483       ncont++;
01484       i = i+spec_bin;
01485     }
01486     else
01487       i++;
01488       }
01489     }
01490     else 
01491       i++;
01492   }
01493 
01494   sinfo_msg ("Number of bins centered on emission features:");
01495   sinfo_msg ("  %i - %i (merged) = %i", norig, nmerged, nline);
01496   sinfo_msg ("  %i continuum bins", ncont);
01497 
01498   s1 = (int*)cpl_calloc(norig-nmerged+ncont, sizeof(int));
01499   s2 = (int*)cpl_calloc(norig-nmerged+ncont, sizeof(int));
01500 
01501 
01502   /* 
01503    *  Merge arrays sorted 
01504    */
01505   i=0;
01506   j=0;
01507   k=0;
01508   while (k<norig-nmerged+ncont) {
01509     if (i<norig && j<ncont && x1b[i]<x1[j]) {
01510       s1[k] = x1b[i];
01511       s2[k] = x2b[i];
01512       k++;
01513       i++;
01514     }
01515     else if (i<norig && j<ncont && x1b[i]>x1[j]) {
01516       s1[k] = x1[j];
01517       s2[k] = x2[j];
01518       k++;
01519       j++;
01520     }
01521     else if (i == norig) {
01522       s1[k] = x1[j];
01523       s2[k] = x2[j];
01524       k++;
01525       j++;
01526     }
01527     else if (j == ncont) {
01528       s1[k] = x1b[i];
01529       s2[k] = x2b[i];
01530       k++;
01531       i++;
01532     }
01533     else {
01534       /* Should never happen */
01535       sinfo_msg_warning("Something went wrong when combining "
01536                         "the bins %i and %i", i,j);
01537       break;
01538     }
01539   }
01540 
01541   for (i=0; i<nline+ncont; i++)
01542     sinfo_msg_debug ("Bin start: %i end %i", s1[i], s2[i]);
01543 
01544   *start = s1;
01545   *end   = s2;
01546 
01547   cpl_free (pos);
01548   cpl_free (x1);
01549   cpl_free (x2);
01550   cpl_free (x1b);
01551   cpl_free (x2b);
01552   cpl_free (flux);
01553   cpl_free (spec);
01554   cpl_free (spec_line);
01555   cpl_free (spec_cont);
01556 
01557   return (nline+ncont);
01558 }
01559 
01560 
01561 /*-------------------------------------------------------------------------*/
01588 /*--------------------------------------------------------------------------*/
01589 
01590 static int
01591 sinfo_juha_function1d_natural_spline(
01592     double  *      x,
01593     double  *      y,
01594     int            len,
01595     double  *      splX,
01596     double  *      splY,
01597     int            splLen
01598 )
01599 {
01600   int                         end;
01601   int                         loc, found;
01602   register int        i, j, n;
01603   double              *       h;      /* vector of deltas in x */
01604   double              *       alpha;
01605   double              *       l,
01606     *       mu,
01607     *       z,
01608     *       a,
01609     *       b,
01610     *       c,
01611     *       d,
01612     v;
01613 
01614   end = len - 1;
01615   
01616   a = cpl_malloc(sizeof(double) * splLen * 9) ;
01617   b = a + len;
01618   c = b + len;
01619   d = c + len;
01620   h = d + len;
01621   l = h + len;
01622   z = l + len;
01623   mu = z + len;
01624   alpha = mu + len;
01625   
01626   for (i = 0; i < len; i++) {
01627     a[i] = (double)y[i];
01628   }
01629   
01630   /* Calculate vector of differences */
01631   for (i = 0; i < end; i++) {
01632     h[i] = (double)x[i + 1] - (double)x[i];
01633     if (h[i] < 0.0) {
01634       cpl_free(a) ;
01635       return -1;
01636     }
01637   }
01638   
01639   /* Calculate alpha vector */
01640   for (n = 0, i = 1; i < end; i++, n++) {
01641     /* n = i - 1 */
01642     alpha[i] = 3.0 * ((a[i+1] / h[i]) - (a[i] / h[n]) - (a[i] / h[i]) +
01643               (a[n] / h[n]));
01644   }
01645   
01646   /* Vectors to solve the tridiagonal matrix */
01647   l[0] = l[end] = 1.0;
01648   mu[0] = mu[end] = 0.0;
01649   z[0] = z[end] = 0.0;
01650   c[0] = c[end] = 0.0;
01651   
01652   /* Calculate the intermediate results */
01653   for (n = 0, i = 1; i < end; i++, n++) {
01654     /* n = i-1 */
01655     l[i] = 2 * (h[i] + h[n]) - h[n] * mu[n];
01656     mu[i] = h[i] / l[i];
01657     z[i] = (alpha[i] - h[n] * z[n]) / l[i];
01658   }
01659   for (n = end, j = end - 1; j >= 0; j--, n--) {
01660     /* n = j + 1 */
01661     c[j] = z[j] - mu[j] * c[n];
01662     b[j] = (a[n] - a[j]) / h[j] - h[j] * (c[n] + 2.0 * c[j]) / 3.0;
01663     d[j] = (c[n] - c[j]) / (3.0 * h[j]);
01664   }
01665   
01666   /* Now calculate the new values */
01667   for (j = 0; j < splLen; j++) {
01668     v = (double)splX[j];
01669     splY[j] = (float)0;
01670     
01671     /* Is it outside the interval? */
01672     if ((v < (double)x[0]) || (v > (double)x[end])) {
01673       continue;
01674     }
01675     /* Search for the interval containing v in the x vector */
01676     loc = sinfo_function1d_search_value(x, len, (double)v, &found);
01677     if (found) {
01678       splY[j] = y[loc];
01679     } else {
01680       loc--;
01681       v -= (double)x[loc];
01682       splY[j] = (float)( a[loc] + v * (b[loc] + v * (c[loc] + v * d[loc])));
01683     }
01684   }
01685   cpl_free(a) ;
01686   return 0;
01687 }
01688 
01689 /*-------------------------------------------------------------------------*/
01705 /*--------------------------------------------------------------------------*/
01706 
01707 static int
01708 sinfo_function1d_search_value(
01709     double  *       x,
01710     int             len,
01711     double          key,
01712     int     *       foundPtr
01713 )
01714 {
01715     int high,
01716                 low,
01717                 middle;
01718 
01719     low  = 0;
01720     high = len - 1;
01721 
01722     while (high >= low) {
01723       middle = (high + low) / 2;
01724       if (key > x[middle]) {
01725     low = middle + 1;
01726       } else if (key < x[middle]) {
01727     high = middle - 1;
01728       } else {
01729     *foundPtr = 1;
01730     return (middle);
01731       }
01732     }
01733     *foundPtr = 0;
01734     return (low);
01735 }
01736 
01737 
01738 /*
01739 cpl_vector * sinfo_vector_filter_median_create(
01740         const cpl_vector * v, 
01741         int                hw)
01742 {
01743     cpl_vector * filtered;
01744     double     * row;
01745     int          i, j, k, size;
01746     double temp;
01747     
01748     size = cpl_vector_get_size(v);
01749     filtered = cpl_vector_new(size);
01750 
01751     row = cpl_malloc((2*hw+1) * sizeof(double));
01752     for (i=0; i<size; i++) {
01753       k = 0;
01754       for (j=-hw; j<=hw; j++) 
01755     if ( (i+j) >= 0 && (i+j) < size) {
01756       temp = cpl_vector_get (v, i+j);
01757       row[k] = temp;
01758       k++;
01759     }
01760       cpl_tools_sort_double (row, k);
01761       if (k%2 == 1)
01762     temp = row[k/2];
01763       else
01764     temp = row[k/2-1];
01765       cpl_vector_set (filtered, i, temp);
01766     }
01767     cpl_free(row);
01768     return filtered;
01769 }
01770 */
01771 
01772 static cpl_vector * 
01773 sinfo_juha_vector_filter_median_create(
01774         const cpl_vector * v, 
01775         int                hw)
01776 {
01777     cpl_vector * filtered=NULL;
01778     double     * row=NULL;
01779     int          i, j, k,  size;
01780     double temp;
01781     
01782     size = cpl_vector_get_size(v);
01783     filtered = cpl_vector_new(size);
01784 
01785     row = cpl_malloc((2*hw+1) * sizeof(double));
01786     for (i=0; i<size; i++) {
01787       k = 0;
01788       for (j=-hw; j<=hw; j++) 
01789     if ( (i+j) >= 0 && (i+j) < size) {
01790       temp = cpl_vector_get (v, i+j);
01791       row[k] = temp;
01792       k++;
01793     }
01794       sinfo_tools_sort_double (row, k);
01795 
01796       if (k%2 == 1)
01797     temp = row[k/2];
01798       else
01799     temp = row[k/2-1];
01800       cpl_vector_set (filtered, i, temp);
01801     }
01802     cpl_free(row);
01803     return filtered;
01804 }
01805 
01806 #define CPL_PIX_STACK_SIZE 50
01807 /*---------------------------------------------------------------------------*/
01818 /*---------------------------------------------------------------------------*/
01819 static cpl_error_code sinfo_tools_sort_double(
01820         double  *   pix_arr,
01821         int         n)
01822 {
01823     int         i, ir, j, k, l;
01824     int     *   i_stack ;
01825     int         j_stack ;
01826     double      a ;
01827     
01828     /* Check entries */
01829     cpl_ensure(pix_arr, CPL_ERROR_NULL_INPUT, CPL_ERROR_NULL_INPUT) ;
01830 
01831     ir = n ;
01832     l = 1 ;
01833     j_stack = 0 ;
01834     i_stack = malloc(CPL_PIX_STACK_SIZE * sizeof(double)) ;
01835     for (;;) {
01836         if (ir-l < 7) {
01837             for (j=l+1 ; j<=ir ; j++) {
01838                 a = pix_arr[j-1];
01839                 for (i=j-1 ; i>=1 ; i--) {
01840                     if (pix_arr[i-1] <= a) break;
01841                     pix_arr[i] = pix_arr[i-1];
01842                 }
01843                 pix_arr[i] = a;
01844             }
01845             if (j_stack == 0) break;
01846             ir = i_stack[j_stack-- -1];
01847             l  = i_stack[j_stack-- -1];
01848         } else {
01849             k = (l+ir) >> 1;
01850             SINFO_DOUBLE_SWAP(pix_arr[k-1], pix_arr[l])
01851             if (pix_arr[l] > pix_arr[ir-1]) {
01852                 SINFO_DOUBLE_SWAP(pix_arr[l], pix_arr[ir-1])
01853             }
01854             if (pix_arr[l-1] > pix_arr[ir-1]) {
01855                 SINFO_DOUBLE_SWAP(pix_arr[l-1], pix_arr[ir-1])
01856             }
01857             if (pix_arr[l] > pix_arr[l-1]) {
01858                 SINFO_DOUBLE_SWAP(pix_arr[l], pix_arr[l-1])
01859             }
01860             i = l+1;
01861             j = ir;
01862             a = pix_arr[l-1];
01863             for (;;) {
01864                 do i++; while (pix_arr[i-1] < a);
01865                 do j--; while (pix_arr[j-1] > a);
01866                 if (j < i) break;
01867                 SINFO_DOUBLE_SWAP(pix_arr[i-1], pix_arr[j-1]);
01868             } 
01869             pix_arr[l-1] = pix_arr[j-1];
01870             pix_arr[j-1] = a; 
01871             j_stack += 2; 
01872             if (j_stack > CPL_PIX_STACK_SIZE) {
01873                 /* Should never reach here */
01874                 free(i_stack);
01875                 return CPL_ERROR_ILLEGAL_INPUT ; 
01876             } 
01877             if (ir-i+1 >= j-l) {
01878                 i_stack[j_stack-1] = ir;
01879                 i_stack[j_stack-2] = i;
01880                 ir = j-1;
01881             } else {
01882                 i_stack[j_stack-1] = j-1;
01883                 i_stack[j_stack-2] = l;
01884                 l = i;
01885             }
01886         }
01887     }
01888     free(i_stack) ;
01889     return CPL_ERROR_NONE ;
01890 }
01891 
01892 static cpl_vector * 
01893 sinfo_vector_filter_median_create(
01894         const cpl_vector * v, 
01895         int                hw)
01896 {
01897     cpl_vector * filtered;
01898     cpl_vector * row;
01899     int          i, j, k, size;
01900     double temp;
01901     
01902     /* Create the filtered vector */
01903     size = cpl_vector_get_size(v);
01904     filtered = cpl_vector_new(size);
01905 
01906     /* median filter on all central items */
01907     row = cpl_vector_new((2*hw+1));
01908     for (i=0; i<size; i++) {
01909       k = 0;
01910       for (j=-hw; j<=hw; j++) 
01911     if ( (i+j) >= 0 && (i+j) < size) {
01912       temp = cpl_vector_get (v, i+j);
01913       cpl_vector_set(row,k,temp);
01914       k++;
01915     }
01916       /*  this returns ~2e8 when all the values are 1.0....*/
01917 /*       temp = cpl_tools_get_median_double(row, k); */
01918       cpl_vector_sort(row, +1);
01919       if (k%2 == 1) {
01920     temp = cpl_vector_get(row,k/2);
01921       }
01922       else {
01923     temp = cpl_vector_get(row,k/2-1);
01924       }
01925 
01926       sinfo_msg("value = %g ", temp);
01927       cpl_vector_set (filtered, i, temp);
01928     }
01929     cpl_vector_delete(row);
01930     return filtered;
01931 }
01932 
01933 /*
01934  *  A NaN safe version of cpl_image_get_median_window
01935  */
01936 static double 
01937 sinfo_image_get_median_window  (const cpl_image *image, 
01938                        int llx, int lly, int urx, int ury)
01939 {
01940   cpl_image *window;
01941   float *data;
01942   double *array, median;
01943   int size, i,j;
01944 
01945   window = cpl_image_extract (image, llx, lly, urx, ury);
01946   size = (urx-llx+1)*(ury-lly+1);
01947   data = cpl_image_get_data_float(window);
01948 
01949   array = (double*)cpl_calloc ( size, sizeof(double));
01950   j = 0;
01951   for (i=0; i<size; i++)
01952     if (!isnan(data[i]))
01953       array[j++] = data[i];
01954   
01955   if (j>0)
01956     sinfo_tools_sort_double (array, j);
01957   
01958   if (j == 0 || 2*j<size)
01959     median = 0./0.;
01960   else if (j%2 == 1) 
01961     median = array[j/2];
01962   else 
01963     median = array[j/2-1];
01964 
01965   cpl_image_delete (window);
01966   cpl_free (array);
01967   
01968   return (median);
01969 }

Generated on Wed Jan 17 08:33:44 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4