irplib_calib.c

00001 /* $Id: irplib_calib.c,v 1.14 2006/11/16 12:26:51 yjung Exp $
00002  *
00003  * This file is part of the irplib package
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program 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
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: yjung $
00023  * $Date: 2006/11/16 12:26:51 $
00024  * $Revision: 1.14 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 /*-----------------------------------------------------------------------------
00033                                    Includes
00034  -----------------------------------------------------------------------------*/
00035 
00036 #include <math.h>
00037 #include <cpl.h>
00038 
00039 #include "irplib_calib.h"
00040 
00041 
00042 static int  
00043 irplib_get_clean_mean_window(cpl_image* img, 
00044                              const int llx, 
00045                              const int lly, 
00046                              const int urx, int ury, 
00047                              const int kappa, 
00048                              const int nclip, 
00049                              double* clean_mean, 
00050                              double* clean_stdev);
00051 
00052 static double irplib_pfits_get_dit(const cpl_propertylist * plist);
00053 static double irplib_pfits_get_exp_time(const cpl_propertylist* plist);
00054 /*----------------------------------------------------------------------------*/
00058 /*----------------------------------------------------------------------------*/
00059 
00061 /*---------------------------------------------------------------------------*/
00068 /*---------------------------------------------------------------------------*/
00069 static double irplib_pfits_get_dit(const cpl_propertylist * plist)
00070 {
00071     return cpl_propertylist_get_double(plist,"ESO DET DIT");
00072 }
00073 
00074 /*---------------------------------------------------------------------------*/
00080 /*---------------------------------------------------------------------------*/
00081 static double irplib_pfits_get_exp_time(const cpl_propertylist* plist)
00082 {
00083   
00084     return cpl_propertylist_get_double(plist,"EXPTIME");
00085 
00086 }
00087 
00088 
00104 static int  
00105 irplib_get_clean_mean_window(cpl_image* img, 
00106                              const int llx, 
00107                              const int lly, 
00108                              const int urx, int ury, 
00109                              const int kappa, 
00110                              const int nclip, 
00111                              double* clean_mean, 
00112                              double* clean_stdev)
00113 {
00114 
00115 
00116   double mean=0;
00117   double stdev=0;
00118   double threshold=0;
00119   double lo_cut=0;
00120   double hi_cut=0;
00121   cpl_mask* mask=NULL;
00122   cpl_image* tmp=NULL;
00123   cpl_stats* stats=NULL;
00124   int i=0;
00125   
00126   tmp=cpl_image_extract(img,llx,lly,urx,ury);
00127   cpl_image_accept_all(tmp);
00128   for(i=0;i<nclip;i++) {
00129    
00130 
00131     cpl_stats_delete(stats);
00132     stats = cpl_stats_new_from_image(tmp, CPL_STATS_MEAN | CPL_STATS_STDEV);
00133     mean = cpl_stats_get_mean(stats);
00134     stdev = cpl_stats_get_stdev(stats);
00135 
00136     threshold=kappa*stdev;
00137     lo_cut=mean-threshold;
00138     hi_cut=mean+threshold;
00139 
00140     cpl_image_accept_all(tmp);
00141     mask=cpl_mask_threshold_image_create(tmp,lo_cut,hi_cut);
00142 
00143     cpl_mask_not(mask);
00144     cpl_image_reject_from_mask(tmp,mask);
00145     cpl_mask_delete(mask);
00146 
00147 
00148   }
00149   *clean_mean=mean;
00150   *clean_stdev=stdev;
00151   cpl_image_delete(tmp);
00152   cpl_stats_delete(stats);
00153  
00154   return 0;
00155 
00156 
00157 }
00158 
00159 /*---------------------------------------------------------------------------*/
00176 /*---------------------------------------------------------------------------*/
00177 
00178 
00179 cpl_table* 
00180 irplib_compute_gain(
00181                 cpl_frameset* son, 
00182                 cpl_frameset* sof, 
00183                 int* zone1, 
00184                 int* zone2,  
00185                 const int kappa,
00186                 const int nclip)
00187 {
00188 
00189   cpl_frame*    frm=NULL;
00190 
00191   cpl_image* img_on1=NULL;
00192   cpl_image* img_on2=NULL;
00193   cpl_image* img_on_dif=NULL;
00194 
00195   cpl_image* img_of1=NULL;
00196   cpl_image* img_of2=NULL;
00197   cpl_image* img_of_dif=NULL;
00198 
00199   cpl_table* res_tbl=NULL;
00200   cpl_vector* dit_on=NULL;
00201   cpl_vector* dit_of=NULL;
00202   cpl_vector* exptime_on=NULL;
00203   cpl_vector* exptime_of=NULL;
00204   cpl_propertylist* plist=NULL;
00205 
00206   int non=0;
00207   int nof=0;
00208   int nfr=0;
00209   int llx;
00210   int lly;
00211   int urx;
00212   int ury;
00213 
00214   double avg_on1=0;
00215   double avg_on2=0;
00216   double avg_of1=0;
00217   double avg_of2=0;
00218   double std=0;
00219 
00220   double sig_on_dif=0;
00221   double sig_of_dif=0;
00222   char* name=NULL;
00223   int i=0;
00224   int m=0;
00225 
00226   double gain=0;
00227   double dit_ref=0;
00228   double dit_tmp=0;
00229   double exptime_ref=0;
00230   double exptime_tmp=0;
00231 
00232  
00233   non = cpl_frameset_get_size(son);
00234   nof = cpl_frameset_get_size(sof);
00235   nfr = (non <= nof) ? non : nof;
00236 
00237   dit_on=cpl_vector_new(nfr);
00238   dit_of=cpl_vector_new(nfr);  
00239   exptime_on=cpl_vector_new(nfr);
00240   exptime_of=cpl_vector_new(nfr);
00241 
00242   for(i=0;i<nfr;i++) {
00243 
00244     frm=cpl_frameset_get_frame(son,i);
00245     name=(char*)cpl_frame_get_filename(frm);
00246     plist=cpl_propertylist_load(name,0);
00247     dit_ref=irplib_pfits_get_dit(plist);
00248     exptime_ref=(double)irplib_pfits_get_exp_time(plist);
00249     cpl_propertylist_delete(plist);
00250     cpl_vector_set(dit_on,i,dit_ref);
00251     cpl_vector_set(exptime_on,i,exptime_ref);
00252 
00253     frm=cpl_frameset_get_frame(sof,i);
00254     name=(char*)cpl_frame_get_filename(frm);
00255     plist=cpl_propertylist_load(name,0);
00256     dit_ref=irplib_pfits_get_dit(plist);
00257     exptime_ref=(double)irplib_pfits_get_exp_time(plist);
00258     cpl_propertylist_delete(plist);
00259     cpl_vector_set(dit_of,i,dit_ref);
00260     cpl_vector_set(exptime_of,i,exptime_ref);
00261 
00262   }
00263 
00264 
00265   llx=zone1[0];
00266   lly=zone1[1];
00267   urx=zone1[2];
00268   ury=zone1[3];
00269 
00270 
00271 
00272   res_tbl=cpl_table_new(nfr);
00273   cpl_table_new_column(res_tbl,"adu", CPL_TYPE_DOUBLE);
00274   cpl_table_new_column(res_tbl,"gain", CPL_TYPE_DOUBLE);
00275  
00276   for(i=0;i<nfr;i++) {
00277     frm=cpl_frameset_get_frame(son,i);
00278     name=(char*)cpl_frame_get_filename(frm);
00279     img_on1=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00280 
00281     frm=cpl_frameset_get_frame(sof,i);
00282     name=(char*)cpl_frame_get_filename(frm);
00283     img_of1=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00284 
00285 
00286     dit_ref=cpl_vector_get(dit_on,i);
00287     exptime_ref=cpl_vector_get(exptime_on,i);
00288 
00289    
00290     for(m=0;m<nfr; m++) {
00291       if(m != i) {
00292     frm=cpl_frameset_get_frame(son,m);
00293     name=(char*)cpl_frame_get_filename(frm);
00294     dit_tmp=cpl_vector_get(dit_on,m);
00295     exptime_tmp=cpl_vector_get(exptime_on,m);
00296     if(dit_tmp == dit_ref && exptime_tmp == exptime_ref) {
00297       img_on2=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00298       frm=cpl_frameset_get_frame(sof,m);
00299       name=(char*)cpl_frame_get_filename(frm);
00300       img_of2=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00301 
00302       img_on_dif=cpl_image_subtract_create(img_on1,img_on2);
00303       img_of_dif=cpl_image_subtract_create(img_of1,img_of2);
00304       
00305       irplib_get_clean_mean_window(img_on1,llx,lly,urx,ury,kappa,
00306                                       nclip,&avg_on1,&std);
00307       irplib_get_clean_mean_window(img_on2,llx,lly,urx,ury,kappa,
00308                                       nclip,&avg_on2,&std);
00309       irplib_get_clean_mean_window(img_of1,llx,lly,urx,ury,kappa,
00310                                       nclip,&avg_of1,&std);
00311       irplib_get_clean_mean_window(img_of2,llx,lly,urx,ury,kappa,
00312                                       nclip,&avg_of2,&std);
00313       
00314           cpl_flux_get_noise_window(img_on_dif,zone2,2,100, &sig_on_dif, NULL);
00315           cpl_flux_get_noise_window(img_of_dif,zone2,2,100, &sig_of_dif, NULL);
00316       
00317       cpl_image_delete(img_on2);
00318       cpl_image_delete(img_of2);
00319       cpl_image_delete(img_on_dif);
00320       cpl_image_delete(img_of_dif);
00321 
00322           gain=((avg_on1+avg_on2)-(avg_of1+avg_of2))/
00323                ((sig_on_dif*sig_on_dif)-(sig_of_dif*sig_of_dif));
00324 
00325           cpl_table_set_double(res_tbl,"gain",m,gain);
00326           cpl_table_set_double(res_tbl,"adu",m,
00327                                ((avg_on1+avg_on2)/2-(avg_of1+avg_of2)/2));
00328 
00329     }
00330       }
00331     }
00332     cpl_image_delete(img_on1);
00333     cpl_image_delete(img_of1);
00334   }
00335   
00336 
00337   cpl_vector_delete(dit_on);
00338   cpl_vector_delete(dit_of);
00339   cpl_vector_delete(exptime_on);
00340   cpl_vector_delete(exptime_of);
00341 
00342   return res_tbl;
00343 
00344 }
00345 
00346 /* --------------------------------------------------------------------------*/
00356 /*---------------------------------------------------------------------------*/
00357 
00358 
00359 cpl_table* irplib_compute_linearity(cpl_frameset* son, cpl_frameset* sof)
00360 {
00361 
00362 
00363   cpl_frame*    frm=NULL;
00364 
00365   int* status=0;
00366   int non=0;
00367   int nof=0;
00368   int nfr=0;
00369   int i=0;
00370   double med_on=0;
00371   double avg_on=0;
00372   double med_of=0;
00373   double avg_of=0;
00374   double med_dit=0;
00375   double avg_dit=0;
00376 
00377   double med=0;
00378   double avg=0;
00379 
00380   char* name=NULL;
00381   cpl_image* img=NULL;
00382   cpl_vector* vec_adl=NULL;
00383   cpl_vector* vec_dit=NULL;
00384   cpl_vector* vec_avg=NULL;
00385   cpl_vector* vec_med=NULL;
00386   cpl_vector* vec_avg_dit=NULL;
00387   cpl_vector* vec_med_dit=NULL;
00388   cpl_propertylist* plist=NULL;
00389 
00390   double dit=0;
00391   cpl_table* lin_tbl=NULL;
00392 
00393  
00394   non = cpl_frameset_get_size(son);
00395   nof = cpl_frameset_get_size(sof);
00396   nfr = (non <= nof) ? non : nof;
00397 
00398   lin_tbl=cpl_table_new(nfr);
00399   cpl_table_new_column(lin_tbl,"med", CPL_TYPE_DOUBLE);
00400   cpl_table_new_column(lin_tbl,"avg", CPL_TYPE_DOUBLE);
00401   cpl_table_new_column(lin_tbl,"med_dit", CPL_TYPE_DOUBLE);
00402   cpl_table_new_column(lin_tbl,"avg_dit", CPL_TYPE_DOUBLE);
00403   cpl_table_new_column(lin_tbl,"dit", CPL_TYPE_DOUBLE);
00404   vec_med=cpl_vector_new(nfr);
00405   vec_avg=cpl_vector_new(nfr);
00406   vec_med_dit=cpl_vector_new(nfr);
00407   vec_avg_dit=cpl_vector_new(nfr);
00408   vec_dit=cpl_vector_new(nfr);
00409   vec_adl=cpl_vector_new(nfr);
00410   for(i=0;i<nfr;i++) {
00411     frm=cpl_frameset_get_frame(son,i);
00412     name=(char*)cpl_frame_get_filename(frm);
00413     img=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00414     med_on=cpl_image_get_median(img);
00415     avg_on=cpl_image_get_mean(img);
00416     cpl_image_delete(img);
00417 
00418     frm=cpl_frameset_get_frame(sof,i);
00419     name=(char*)cpl_frame_get_filename(frm);
00420     img=cpl_image_load(name,CPL_TYPE_FLOAT,0,0);
00421     med_of=cpl_image_get_median(img);
00422     avg_of=cpl_image_get_mean(img);
00423     cpl_image_delete(img);
00424     med=med_on-med_of;
00425     avg=avg_on-avg_of;
00426     plist=cpl_propertylist_load(name,0);
00427     dit=(double)irplib_pfits_get_dit(plist);
00428     cpl_propertylist_delete(plist);
00429     avg_dit=avg/dit;
00430     med_dit=med/dit;
00431 
00432     cpl_vector_set(vec_dit,i,dit);
00433     cpl_vector_set(vec_avg,i,avg);
00434     cpl_vector_set(vec_med,i,med);
00435     cpl_vector_set(vec_avg_dit,i,avg_dit);
00436     cpl_vector_set(vec_med_dit,i,med_dit);
00437 
00438 
00439     cpl_table_set_double(lin_tbl,"dit",i,dit);
00440     cpl_table_set_double(lin_tbl,"med",i,med);
00441     cpl_table_set_double(lin_tbl,"avg",i,avg);
00442     cpl_table_set_double(lin_tbl,"med_dit",i,med_dit);
00443     cpl_table_set_double(lin_tbl,"avg_dit",i,avg_dit);
00444 
00445   }
00446   cpl_table_new_column(lin_tbl,"adl", CPL_TYPE_DOUBLE);
00447   med_dit=cpl_vector_get_mean(vec_med_dit);
00448   avg_dit=cpl_vector_get_mean(vec_avg_dit);
00449 
00450   for(i=0;i<nfr;i++) {
00451     dit = cpl_table_get_double(lin_tbl,"dit",i,status);
00452     cpl_vector_set(vec_adl,i,dit*med_dit);
00453     cpl_table_set_double(lin_tbl,"adl",i,dit*med_dit);
00454   }
00455  
00456   
00457   cpl_vector_delete(vec_dit);
00458   cpl_vector_delete(vec_adl);
00459   cpl_vector_delete(vec_avg);
00460   cpl_vector_delete(vec_med);
00461   cpl_vector_delete(vec_avg_dit);
00462   cpl_vector_delete(vec_med_dit);
00463 
00464 
00465   return lin_tbl;
00466 
00467 }
00468 
00469 
00470 /*----------------------------------------------------------------------------*/
00479 /*----------------------------------------------------------------------------*/
00480 int irplib_detlin_correct(
00481         cpl_imagelist       *   ilist,
00482         const char          *   detlin_a,
00483         const char          *   detlin_b,
00484         const char          *   detlin_c)
00485 {
00486     cpl_image       *   ima ;
00487     cpl_image       *   imb ;
00488     cpl_image       *   imc ;
00489     float           *   pima ;
00490     float           *   pimb ;
00491     float           *   pimc ;
00492     float           *   pdata ;
00493     int                 nx, ny, ni ; 
00494     double              coeff_1, coeff_2, val ;
00495     int                 i, j ;
00496 
00497     /* Test entries */
00498     if (!ilist || !detlin_a || !detlin_b || !detlin_c) return -1 ;
00499     
00500     /* Load the 3 coeffs images */
00501     ima = cpl_image_load(detlin_a, CPL_TYPE_FLOAT, 0, 0) ;
00502     imb = cpl_image_load(detlin_b, CPL_TYPE_FLOAT, 0, 0) ;
00503     imc = cpl_image_load(detlin_c, CPL_TYPE_FLOAT, 0, 0) ;
00504     if (!ima || !imb || !imc) {
00505         cpl_msg_error(cpl_func, "Cannot load the detlin images") ;
00506         if (ima) cpl_image_delete(ima) ;
00507         if (imb) cpl_image_delete(imb) ;
00508         if (imc) cpl_image_delete(imc) ;
00509         return -1 ;
00510     }
00511     pima = cpl_image_get_data_float(ima) ;
00512     pimb = cpl_image_get_data_float(imb) ;
00513     pimc = cpl_image_get_data_float(imc) ;
00514     
00515     /* Test sizes */
00516     nx = cpl_image_get_size_x(cpl_imagelist_get(ilist, 0)) ;
00517     ny = cpl_image_get_size_y(cpl_imagelist_get(ilist, 0)) ;
00518     ni = cpl_imagelist_get_size(ilist) ;
00519     if ((cpl_image_get_size_x(ima) != nx) ||
00520             (cpl_image_get_size_x(imb) != nx) ||
00521             (cpl_image_get_size_x(imc) != nx) ||
00522             (cpl_image_get_size_y(ima) != ny) ||
00523             (cpl_image_get_size_y(imb) != ny) ||
00524             (cpl_image_get_size_y(imc) != ny)) {
00525         cpl_msg_error(cpl_func, "Incompatible sizes") ;
00526         cpl_image_delete(ima) ;
00527         cpl_image_delete(imb) ;
00528         cpl_image_delete(imc) ;
00529         return -1 ;
00530     }
00531     
00532     /* Loop on pixels */
00533     for (i=0 ; i<nx*ny ; i++) {
00534         /* Compute the coefficients */
00535         if (fabs(pima[i]) < 1e-30) {
00536             coeff_1 = coeff_2 = (double)0.0 ;
00537         } else {
00538             coeff_1 = (double)pimb[i] / (double)pima[i] ;
00539             coeff_2 = (double)pimc[i] / (double)pima[i] ;
00540         }
00541         /* Correct this pixel in each plane */
00542         for (j=0 ; j<ni ; j++) {
00543             pdata = cpl_image_get_data_float(cpl_imagelist_get(ilist, j)) ;
00544             val = (double)pdata[i] ;
00545             pdata[i]=(float)(val+coeff_1*val*val+coeff_2*val*val*val) ;
00546         }
00547     }
00548     /* Free and return */
00549     cpl_image_delete(ima) ;
00550     cpl_image_delete(imb) ;
00551     cpl_image_delete(imc) ;
00552     return 0 ;
00553 }
00554  
00555 /*----------------------------------------------------------------------------*/
00564 /*----------------------------------------------------------------------------*/
00565 int irplib_flat_dark_bpm_calib(
00566         cpl_imagelist       *   ilist,
00567         const char          *   flat,
00568         const char          *   dark,
00569         const char          *   bpm)
00570 {
00571     cpl_image       *   dark_image ;
00572     cpl_image       *   flat_image ;
00573     cpl_mask        *   bpm_im_bin ;
00574     cpl_image       *   bpm_im_int ;
00575     int                 i ;
00576 
00577     /* Test entries */
00578     if (ilist == NULL) return -1 ;
00579 
00580     /* Dark correction */
00581     if (dark != NULL) {
00582         cpl_msg_info(cpl_func, "Subtract the dark to the images") ;
00583         /* Load the dark image */
00584         if ((dark_image = cpl_image_load(dark, CPL_TYPE_FLOAT, 0, 0)) == NULL) {
00585             cpl_msg_error(cpl_func, "Cannot load the dark %s", dark) ;
00586             return -1 ;
00587         }
00588         /* Apply the dark correction to the images */
00589         if (cpl_imagelist_subtract_image(ilist, dark_image)!=CPL_ERROR_NONE) {
00590             cpl_msg_error(cpl_func, "Cannot apply the dark to the images") ;
00591             cpl_image_delete(dark_image) ;
00592             return -1 ;
00593         }
00594         cpl_image_delete(dark_image) ;
00595     }
00596 
00597     /* Flat-field correction */
00598     if (flat != NULL) {
00599         cpl_msg_info(cpl_func, "Divide the images by the flatfield") ;
00600         /* Load the flat image */
00601         if ((flat_image = cpl_image_load(flat, CPL_TYPE_FLOAT, 0, 0)) == NULL) {
00602             cpl_msg_error(cpl_func, "Cannot load the flat field %s", flat) ;
00603             return -1 ;
00604         }
00605         /* Apply the flatfield correction to the images */
00606         if (cpl_imagelist_divide_image(ilist, flat_image)!=CPL_ERROR_NONE) {
00607             cpl_msg_error(cpl_func, "Cannot apply the flatfield to the images") ;
00608             cpl_image_delete(flat_image) ;
00609             return -1 ;
00610         }
00611         cpl_image_delete(flat_image) ;
00612     }
00613 
00614     /* Correct the bad pixels if requested */
00615     if (bpm != NULL) {
00616         cpl_msg_info(cpl_func, "Correct the bad pixels in the images") ;
00617          /* Load the bad pixels image */
00618         if ((bpm_im_int = cpl_image_load(bpm, CPL_TYPE_INT, 0, 0)) == NULL) {
00619             cpl_msg_error(cpl_func, "Cannot load the bad pixel map %s", bpm) ;
00620             return -1 ;
00621         }
00622         /* Convert the map from integer to binary */
00623         bpm_im_bin = cpl_mask_threshold_image_create(bpm_im_int, -0.5, 0.5) ;
00624         cpl_mask_not(bpm_im_bin) ;
00625         cpl_image_delete(bpm_im_int) ;
00626         /* Apply the bad pixels cleaning */
00627         for (i=0 ; i<cpl_imagelist_get_size(ilist) ; i++) {
00628             cpl_image_reject_from_mask(cpl_imagelist_get(ilist, i), bpm_im_bin);
00629             if (cpl_detector_interpolate_rejected(
00630                         cpl_imagelist_get(ilist, i)) != CPL_ERROR_NONE) {
00631                 cpl_msg_error(cpl_func, "Cannot clean the bad pixels in obj %d",
00632                         i+1);
00633                 cpl_mask_delete(bpm_im_bin) ;
00634                 return -1 ;
00635             }
00636         }
00637         cpl_mask_delete(bpm_im_bin) ;
00638     }
00639 
00640     /* Return */
00641     return 0 ;
00642 }
00643 

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