sinfo_new_nst.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*----------------------------------------------------------------------------
00020  
00021      File name    :     spiffi_north_south_test.c
00022    Author       :   A. Modigliani
00023    Created on   :   Sep 17, 2003
00024    Description  : 
00025 
00026  * nsh.h
00027  * Result of a north-south test exposure are 32 continuum spectra of a 
00028  * pinhole that means one spectrum in each slitlet at the same spatial 
00029  * position.
00030  * Each spectrum is fitted in sp[atial direction by a Gaussian to get the 
00031  * sub-pixel positions for each row.
00032  *
00033  * Then the distances are determined in each row and averaged
00034  *
00035  * Result: are distances of each slitlet from each other => 31 values stored 
00036  *  in an ASCII file this Python script needs a frame of a pinhole source with 
00037  *  a continuous spectrum, that is shifted exactly perpendicular to the 
00038  *  slitlets. It fits the spectra in spatial direction by a Gaussian fit 
00039  *  function and therefore determines the sub-pixel position of the source.
00040  *
00041  *  Then the distances of the slitlets from each other are determined and 
00042  *   saved in an ASCII list. 
00043  *
00044 
00045 
00046  ---------------------------------------------------------------------------*/
00047 
00048 #ifdef HAVE_CONFIG_H
00049 #  include <config.h>
00050 #endif
00051 /*----------------------------------------------------------------------------
00052                                 Includes
00053  ---------------------------------------------------------------------------*/
00054 #include "sinfo_new_nst.h"
00055 #include "sinfo_pro_save.h"
00056 #include "sinfo_pro_types.h"
00057 #include "sinfo_functions.h"
00058 #include "sinfo_ns_ini_by_cpl.h"
00059 #include "sinfo_cube_construct.h"
00060 #include "sinfo_utilities.h"
00061 #include "sinfo_utils_wrappers.h"
00062 #include "sinfo_error.h"
00063 #include "sinfo_globals.h"
00064 
00065 /*----------------------------------------------------------------------------
00066                                 Defines
00067  ---------------------------------------------------------------------------*/
00068 
00076 /*----------------------------------------------------------------------------
00077                              Function Definitions
00078  ---------------------------------------------------------------------------*/
00079 
00080 /*----------------------------------------------------------------------------
00081    Function     :       north_south_test()
00082    In           :       ini_file: file name of according .ini file
00083    Out          :       integer (0 if it worked, -1 if it doesn't) 
00084    Job          : Result of a north-south test exposure are 32 continuum 
00085                   spectra of a pinhole  that means one spectrum in each 
00086                   slitlet at the same spatial position.
00087 
00088                   Each spectrum is fitted in sp[atial direction by a Gaussian
00089                   to get the sub-pixel positions for each row.
00090 
00091                   Then the distances are determined in each row and averaged
00092 
00093                   Result: are distances of each slitlet from each other =>
00094                           31 values stored in an ASCII file
00095 
00096 
00097  ---------------------------------------------------------------------------*/
00098 int 
00099 sinfo_new_nst(const char* plugin_id, 
00100               cpl_parameterlist* config, 
00101               cpl_frameset* sof)
00102 {
00103   ns_config * cfg=NULL ;
00104   cpl_imagelist * list_object=NULL ;
00105   cpl_imagelist * list_off=NULL;
00106   cpl_image * im_on=NULL;
00107   cpl_image * im_on_sub=NULL ;
00108   cpl_image * im_on_ind=NULL ;
00109   cpl_image * im_on_gauss=NULL ;
00110   cpl_image * im_mask=NULL ;
00111   cpl_image * im_off=NULL ;
00112 
00113   char* name=NULL;
00114   char tbl_name[MAX_NAME_SIZE];
00115   int typ=0;
00116   int i = 0;
00117 
00118   int nob =0;
00119   int nof =0;
00120 
00121   float* distances=NULL;
00122   cpl_table* tbl_dist=NULL;
00123 
00124   cpl_vector* qc_dist=NULL;
00125 
00126   double  qc_dist_mean=0;
00127   double  qc_dist_stdev=0;
00128   cpl_frameset* raw=NULL;
00129 
00130   cpl_table* qclog_tbl=NULL;
00131   char key_value[MAX_NAME_SIZE];
00132   char key_name[MAX_NAME_SIZE];
00133   int  no=0;
00134   double lo_cut=0.;
00135   double hi_cut=0.;
00136 
00137   /* 
00138        -----------------------------------------------------------------
00139        1) parse the file names and parameters to the ns_config data 
00140           structure cfg
00141        -----------------------------------------------------------------
00142   */
00143 
00144   /* 
00145       parse the file names and parameters to the ns_config data structure cfg 
00146   */
00147 
00148   sinfo_msg("Parse cpl input");
00149   check_nomsg(raw=cpl_frameset_new());
00150   cknull(cfg = sinfo_parse_cpl_input_ns(config,sof,&raw),
00151      "could not parse cpl input!") ;
00152 
00153   if (cfg->maskInd == 1) {
00154     if(sinfo_is_fits_file(cfg->mask) != 1) {
00155       sinfo_msg_error("Input file %s is not FITS",cfg->mask); 
00156       goto cleanup;
00157     }
00158   }
00159   /*
00160      --------------------------------------------------------------------
00161      stack the frames in data cubes and take the clean mean of all frames
00162      --------------------------------------------------------------------
00163   */
00164   /* allocate memory for lists of object and off-frames */
00165   sinfo_msg("stack the frames in data cubes");
00166   check(list_object = cpl_imagelist_new(),
00167           "could not allocate memory");
00168 
00169   if (cfg->noff > 0 ) {
00170     check(list_off = cpl_imagelist_new(),
00171       "could not allocate memory");
00172   }
00173  
00174   /*
00175       #build different image lists for the different cases----
00176   */
00177   sinfo_msg("build different image lists for the different cases");
00178   nob = 0;
00179   nof = 0;
00180 
00181   for (i=0; i< cfg->nframes; i++){
00182     name = cfg->framelist[i];
00183     if(sinfo_is_fits_file(name) != 1) {
00184       sinfo_msg_error("Input file %s is not FITS",name);
00185       goto cleanup;
00186     } else {
00187       typ = sinfo_new_intarray_get_value( cfg->frametype, i );
00188       if (typ == 1) {
00189     cpl_imagelist_set(list_object,
00190                           cpl_image_load(name,CPL_TYPE_FLOAT,0,0),nob);
00191     nob = nob + 1;
00192       } else {
00193     cpl_imagelist_set(list_off,
00194                           cpl_image_load(name,CPL_TYPE_FLOAT,0,0),nof);
00195     nof = nof + 1;
00196       }
00197     }
00198   }
00199 
00200   if (cfg->noff != nof || cfg->nobj != nob ){       
00201       sinfo_msg_error("something wrong with the number of the "
00202                       "different types of frames");
00203       goto cleanup;
00204   }
00205   
00206   /*
00207   #---take the average of the different cubes -------------
00208   */
00209   sinfo_msg("take the average of the different cubes");
00210 
00211   check_nomsg(no=cpl_imagelist_get_size(list_object));
00212   lo_cut=(floor)(cfg->loReject*no+0.5);
00213   hi_cut=(floor)(cfg->hiReject*no+0.5);
00214   check(im_on=cpl_imagelist_collapse_minmax_create(list_object,lo_cut,hi_cut),
00215               "sinfo_average_with_rejection failed" );
00216 
00217   if (cfg->noff != 0) {
00218     /*
00219       im_off = sinfo_average_with_rejection( cube_off, 
00220                                        cfg->loReject, cfg->hiReject );
00221     */
00222     check_nomsg(no=cpl_imagelist_get_size(list_off));
00223     lo_cut=(floor)(cfg->loReject*no+0.5);
00224     hi_cut=(floor)(cfg->hiReject*no+0.5);
00225     check(im_off=cpl_imagelist_collapse_minmax_create(list_off,lo_cut,hi_cut),
00226       "sinfo_average_with_rejection failed" );
00227       sinfo_free_imagelist(&list_off);
00228   }
00229   sinfo_free_imagelist(&list_object);
00230 
00231   /*
00232   #finally, subtract off from on frames and store the result in the object cube
00233   */
00234 
00235   if (cfg->noff != 0) {
00236     sinfo_msg("subtract off from on frames");
00237     check(im_on_sub = cpl_image_subtract_create(im_on, im_off),
00238       "sinfo_sub_image failed" );
00239 
00240     sinfo_free_image(&im_on);
00241     sinfo_free_image(&im_off);
00242   } else {
00243     check_nomsg(im_on_sub = cpl_image_duplicate(im_on));
00244     sinfo_free_image(&im_on);
00245   }
00246 
00247   /*
00248   #---------------------------------------------------------
00249   # convolution with Gaussian if recommended
00250   #---------------------------------------------------------
00251   */
00252   if (cfg->gaussInd == 1) {
00253     sinfo_msg("convolution with Gaussian");
00254     cknull(im_on_gauss = sinfo_new_convolve_ns_image_by_gauss(im_on_sub, 
00255                                                               cfg->hw),
00256                      "could not carry out sinfo_convolveNSImageByGauss" );
00257 
00258     sinfo_free_image(&im_on_sub);
00259     check_nomsg(im_on_sub = cpl_image_duplicate(im_on_gauss));
00260     sinfo_free_image(&im_on_gauss);
00261   }
00262 
00263   /*
00264   #---------------------------------------------------------
00265   # static bad pixel indication
00266   #---------------------------------------------------------
00267   */
00268 
00269    if (cfg->maskInd == 1) {
00270      sinfo_msg("static bad pixel indication");
00271      check(im_mask = cpl_image_load(cfg->mask,CPL_TYPE_FLOAT,0,0),
00272        "could not load static bad pixel mask" );
00273      cknull(im_on_ind = sinfo_new_mult_image_by_mask(im_on_sub, im_mask),
00274         "could not carry out sinfo_multImageByMask" );
00275      sinfo_free_image(&im_mask);
00276      sinfo_free_image(&im_on_sub);
00277   
00278    } else {
00279       check_nomsg(im_on_ind = cpl_image_duplicate(im_on_sub));
00280       sinfo_free_image(&im_on_sub);
00281    }
00282 
00283    ck0(sinfo_pro_save_ima(im_on_ind,raw,sof,cfg->fitsname,
00284               PRO_SLIT,NULL,plugin_id,config),
00285        "cannot save ima %s", cfg->fitsname);
00286  
00287    /*
00288    #---------------------------------------------------------
00289    # do the north - south - test
00290    #---------------------------------------------------------
00291    */
00292    sinfo_msg("Do the north - south - test");
00293 
00294    cknull(distances = sinfo_north_south_test(im_on_ind, 
00295                        cfg->nslits,
00296                        cfg->halfWidth,
00297                        cfg->fwhm ,
00298                        cfg->minDiff,
00299                        cfg->estimated_dist,
00300                                            cfg->devtol, 
00301                                            IMA_PIX_START, 
00302                                            IMA_PIX_END ),
00303       "North South Test distance determination failed");
00304    sinfo_free_image(&im_on_ind);
00305     /*
00306    sinfo_new_parameter_to_ascii(distances, cfg->nslits - 1, cfg->outName);
00307    */
00308    check_nomsg(tbl_dist = cpl_table_new(cfg->nslits - 1));
00309    check_nomsg(cpl_table_new_column(tbl_dist,"slitlet_distance",CPL_TYPE_FLOAT));
00310    check_nomsg(cpl_table_copy_data_float(tbl_dist,"slitlet_distance",distances));
00311 
00312    strcpy(tbl_name,cfg->outName);
00313 
00314    check_nomsg(qclog_tbl = cpl_table_new(cfg->nslits + 1));
00315    check_nomsg(cpl_table_new_column(qclog_tbl,"key_name", CPL_TYPE_STRING));
00316    check_nomsg(cpl_table_new_column(qclog_tbl,"key_type", CPL_TYPE_STRING));
00317    check_nomsg(cpl_table_new_column(qclog_tbl,"key_value", CPL_TYPE_STRING));
00318    check_nomsg(cpl_table_new_column(qclog_tbl,"key_help", CPL_TYPE_STRING));
00319 
00320    check_nomsg(qc_dist=cpl_vector_new(cfg->nslits - 1));
00321  
00322    for(i=0;i<cfg->nslits - 1;i++) {
00323          snprintf(key_name,MAX_NAME_SIZE-1,"%s%i","QC SL DIST",i);
00324          cpl_table_set_string(qclog_tbl,"key_name",i,key_name);
00325          cpl_table_set_string(qclog_tbl,"key_type",i,"CPL_TYPE_DOUBLE");
00326          snprintf(key_value,MAX_NAME_SIZE-1,"%g",distances[i]);
00327          cpl_table_set_string(qclog_tbl,"key_value",i,key_value);
00328          cpl_table_set_string(qclog_tbl,"key_help",i,"Slitlet distance");
00329 
00330          cpl_vector_set(qc_dist,i,distances[i]);
00331    }
00332    check_nomsg(qc_dist_mean=cpl_vector_get_mean(qc_dist));
00333    check_nomsg(qc_dist_stdev=cpl_vector_get_stdev(qc_dist));
00334 
00335    cpl_table_set_string(qclog_tbl,"key_name",cfg->nslits-1,"QC SL DISTAVG");
00336    cpl_table_set_string(qclog_tbl,"key_type",cfg->nslits-1,"CPL_TYPE_DOUBLE");
00337    snprintf(key_value,MAX_NAME_SIZE-1,"%g",cpl_vector_get_mean(qc_dist));
00338    cpl_table_set_string(qclog_tbl,"key_value",cfg->nslits-1,key_value);
00339    cpl_table_set_string(qclog_tbl,"key_help",cfg->nslits-1,
00340                                   "Average Slitlet distance");
00341 
00342    cpl_table_set_string(qclog_tbl,"key_name",cfg->nslits,"QC SL DISTRMS");
00343    cpl_table_set_string(qclog_tbl,"key_type",cfg->nslits,"CPL_TYPE_DOUBLE");
00344    snprintf(key_value,MAX_NAME_SIZE-1,"%g",qc_dist_stdev);
00345    cpl_table_set_string(qclog_tbl,"key_value",cfg->nslits,key_value);
00346    cpl_table_set_string(qclog_tbl,"key_help",cfg->nslits,
00347                         "RMS Slitlet distance");
00348 
00349    ck0(sinfo_pro_save_tbl(tbl_dist,raw,sof,tbl_name,
00350            PRO_SLITLETS_DISTANCE,qclog_tbl,plugin_id,config),
00351        "cannot dump tbl %s", tbl_name);
00352 
00353    sinfoni_free_vector(&qc_dist);
00354    sinfo_free_table(&tbl_dist);
00355    sinfo_free_table(&qclog_tbl);
00356    sinfo_free_float(&distances);
00357    sinfo_ns_free (&cfg);
00358    sinfo_free_frameset(&raw);
00359 
00360   return 0;
00361 
00362   cleanup:
00363   sinfoni_free_vector(&qc_dist);
00364   sinfo_free_table(&tbl_dist);
00365   sinfo_free_table(&qclog_tbl);
00366   sinfo_free_float(&distances);
00367   sinfo_free_table(&tbl_dist);
00368   sinfo_free_table(&qclog_tbl);
00369   sinfo_free_imagelist(&list_object);
00370   sinfo_free_imagelist(&list_off);
00371   sinfo_free_image(&im_on);
00372   sinfo_free_image(&im_mask);
00373   sinfo_free_image(&im_on_gauss);
00374   sinfo_free_image(&im_on_sub);
00375   sinfo_free_image(&im_off);
00376   sinfo_free_image(&im_on_ind);
00377   sinfo_ns_free (&cfg);
00378   sinfo_free_frameset(&raw);
00379   return -1;
00380 
00381 
00382 }
00383 

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