si_utl_illumcorr.c

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

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