bp_sky.c

00001 /*----------------------------------------------------------------------------
00002    
00003    File name    :       bp_sky.c
00004    Author       :   J. Schreiber
00005    Created on   :   May 5, 2003
00006    Description  :   Generates a bad pixel map associated to a sky frame
00007  ---------------------------------------------------------------------------*/
00008 
00009 /*----------------------------------------------------------------------------
00010                                 Includes
00011  ---------------------------------------------------------------------------*/
00012 #include "bp_sky.h"
00013 #include "badsky_ini_by_cpl.h"
00014 #include "sinfoni_pro_save.h"
00015 
00016 /*----------------------------------------------------------------------------
00017                                 Defines
00018  ---------------------------------------------------------------------------*/
00019 
00020 
00021 /*----------------------------------------------------------------------------
00022                              Function Definitions
00023  ---------------------------------------------------------------------------*/
00024 
00025 /*----------------------------------------------------------------------------
00026    Function     :       bp_sky()
00027    In           :       ini_file: file name of according .ini file
00028    Out          :       integer (0 if it worked, -1 if it doesn't) 
00029    Job          :       
00030 
00031    this function searches for static bad pixels in stacks of flatfield frames, 
00032    that means it stacks the flatfield in a data cube and takes the median 
00033    along the z-axis to be sure to have no cosmics, then the intensity tilt of 
00034    each column is removed. 
00035 
00036    The standard deviation and mean of the pixel values within a defined 
00037    rectangular zone is determined in a way that the extreme low and high 
00038    values are cut off.
00039 
00040    The next step is to indicate each pixel as bad that has a deviation from 
00041    the mean greater than a user defined factor times the standard deviation. 
00042    This step can be leaped over if wished.
00043    The next step is to calculate the median of the 8 nearest neighbors for 
00044    each pixel and to determine the deviation of each pixel value from this 
00045    median. 
00046 
00047    If this deviation is greater than a defined factor times the standard 
00048    deviation the pixel value is replaced by the median. The last step is 
00049    repeated to be able to consider also clusters of bad pixels. Then the 
00050    resulting image is compared with the input image with the column intensity 
00051    tilt removed and each changed pixel is indicated as bad. Finally, a bad 
00052    pixel mask is produced that means each good pixel is marked with 1 and each 
00053    bad pixel with 0.
00054                       
00055  ---------------------------------------------------------------------------*/
00056 
00057 
00058 int bp_sky (cpl_parameterlist* config, 
00059                      cpl_frameset* sof, char* procatg)
00060 {
00061   const char* _id = "bp_sky";
00062     badsky_config * cfg =NULL;
00063 
00064     OneImage ** imagelist;
00065     OneImage ** med ;
00066 
00067     OneImage * medImage =NULL;
00068     OneImage * ima_dark =NULL;
00069     OneImage * medIm =NULL;
00070     OneImage * colImage =NULL;
00071     OneImage * compImage =NULL;
00072     OneImage * maskImage =NULL;
00073     OneImage * threshIm =NULL;
00074 
00075     OneCube * cubeflat =NULL;
00076 
00077     Stats * stats =NULL;
00078 
00079     cpl_image* bp_img=NULL;
00080     /*
00081     cpl_image* ict_img=NULL;
00082     */
00083     cpl_parameter *p=NULL;             
00084 
00085     int i=0;
00086     int n=0;
00087     int half_box_size=0 ;
00088 
00089     cpl_frameset* raw=NULL;
00090     cpl_table* qclog_tbl=NULL;
00091     char* key_value=NULL;
00092 
00093     /* parse the file names and parameters to the bad_config data structure 
00094        cfg */
00095 
00096  
00097     cfg = parse_cpl_input_badsky(config,sof,&raw) ;
00098  
00099     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00100       cpl_msg_error(_id,"error parsing cpl input");
00101       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00102       return -1;
00103     }
00104     if (cfg == NULL)
00105     {
00106         cpl_msg_error (_id,"could not parse cpl input!\n") ;
00107         return -1 ;
00108     }
00109 
00110     /* take a clean mean of the frames */
00111     cpl_msg_info(_id,"take a clean mean of the frames");
00112     imagelist = (OneImage**) cpl_calloc (cfg -> nframes, sizeof(OneImage*)) ;
00113     for ( i = 0 ; i < cfg->nframes ; i++ )
00114     {   
00115         if(is_fits_file (cfg->framelist[i]) != 1) {
00116           cpl_msg_error(_id,"Input file %s is not FITS",cfg->framelist[i] );
00117           return -1;
00118         }
00119         imagelist[i] = load_image(cfg->framelist[i]) ;
00120     }
00121 
00122     cubeflat = list_make_cube(imagelist, cfg-> nframes) ;
00123     if ( cubeflat == NULL )
00124     {
00125         cpl_msg_error(_id," could not make a cube list!\n") ;
00126     return -1 ;
00127     }
00128     
00129     /* finally take the average image of the cube by rejecting the extreme 
00130        values */
00131     medImage = average_with_rejection (cubeflat, cfg->loReject, cfg->hiReject) ;
00132               
00133     if (medImage == NULL)
00134     {
00135         cpl_msg_error(_id," error in average_with_rejection\n") ;
00136         return -1 ;
00137     }
00138 
00139     if (strlen(cfg->dark) > 0) {
00140       cpl_msg_info(_id,"Subtract dark");
00141       ima_dark=load_image(cfg->dark);
00142       medImage = sub_image(medImage,ima_dark);
00143     } else {
00144       cpl_msg_info(_id,"No dark provided");
00145     }
00146     
00147 
00148     /* free memory */
00149     destroy_cube(cubeflat) ;
00150     cpl_free(imagelist) ;
00151 
00152     /*----------------------------------------------
00153      * remove the intensity tilt from every column
00154      * and compute the standard deviation on a rectangular zone
00155      */
00156     medIm = threshImage(medImage, cfg->mincut, cfg->maxcut);
00157     colImage  = colTilt( medIm, cfg->sigmaFactor ) ;
00158     if (colImage == NULL)
00159     {
00160         cpl_msg_error( _id,"colTilt failed" ) ;
00161         return -1 ;
00162     }
00163     stats = imageStatsOnRectangle(colImage, 
00164                                   cfg->loReject, 
00165                                   cfg->hiReject, 
00166                                   cfg->llx, 
00167                                   cfg->lly, 
00168                                   cfg->urx, 
00169                                   cfg->ury) ;
00170     if (stats == NULL)
00171     {
00172         cpl_msg_error (_id," get_image_stats_on_vig failed") ;
00173         return -1 ;
00174     }
00175    
00176     cpl_msg_info (_id, "clean stdev: %f\n", stats->cleanstdev ) ;
00177     cpl_msg_info (_id, "clean mean: %f\n",  stats->cleanmean ) ;
00178     
00179 
00180 
00181     /*
00182     ict_img=cpl_image_wrap_float(colImage->lx,
00183                                 colImage->ly,
00184                                 colImage->data);
00185 
00186 
00187     qclog_tbl = sinfoni_qclog_init(2);
00188     key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00189     sprintf(key_value, "%g",stats->cleanstdev ); 
00190     sinfoni_qclog_add(qclog_tbl,0,(char*) "QC ICTILT STDEV",(char*) "CPL_TYPE_DOUBLE",
00191                       key_value,(char*) "Intensity column clean stdev");
00192 
00193     sprintf(key_value, "%g",stats->cleanmean ); 
00194     sinfoni_qclog_add(qclog_tbl,1,(char*) "QC ICTILT MEAN",(char*) "CPL_TYPE_DOUBLE",
00195                       key_value,(char*) "Intensity column clean mean");
00196 
00197 
00198 
00199     if(-1 == sinfoni_pro_save_ima(ict_img,raw,sof,
00200                (char*)BP_NORM_INT_COL_TILT_CORR_OUT_FILENAME,
00201          (char*) PRO_INT_COL_TILT_COR,qclog_tbl,_id,config)) {
00202          cpl_msg_error(_id,"cannot save ima %s", 
00203            BP_NORM_INT_COL_TILT_CORR_OUT_FILENAME);
00204          cpl_free(key_value);
00205          destroy_image(medImage) ;
00206          destroy_image(medIm) ;
00207          destroy_image(threshIm) ;
00208          if (cfg->threshInd == 1) destroy_image(colImage) ;
00209          if (stats) cpl_free(stats) ;
00210          badsky_free(cfg) ;
00211     
00212          return -1;
00213     }
00214     cpl_table_delete(qclog_tbl);
00215     cpl_free(key_value);
00216     cpl_image_delete(ict_img);
00217     */
00218 
00219     /* indicate pixels with great deviations from the clean mean as bad */
00220     if (cfg->threshInd == 1)
00221     {
00222         threshIm = threshImage(colImage, 
00223                                stats->cleanmean-cfg->meanfactor*stats->cleanstdev, 
00224                                stats->cleanmean+cfg->meanfactor*stats->cleanstdev) ;
00225         if (threshIm == NULL)
00226         {
00227         cpl_msg_error ( _id," threshImage failed" ) ;
00228             return -1 ;   
00229         }
00230     }
00231 
00232     if (cfg->threshInd == 0 )
00233     {
00234         threshIm = colImage ;
00235     }
00236     med = (OneImage**) cpl_calloc (cfg -> iterations, sizeof(OneImage*)) ;
00237     
00238     /* filter iteratively the images by a median filter of the nearest 
00239        neighbors under the condition of a deviation greater than a factor 
00240        times the standard deviation */
00241     cpl_msg_info(_id,"Apply median filter on pixel nearest neighbors");
00242     if (cfg->methodInd == 1) 
00243     {
00244         if (cfg->factor>0)
00245         {
00246             med[0] = medianImage(threshIm, -cfg->factor*stats->cleanstdev) ;
00247             if ( med[0] == NULL )
00248             { 
00249                 cpl_msg_error ( _id," medianImage failed (1)" ) ;
00250                 return -1 ;
00251             }  
00252 
00253             for ( i = 0 ; i < cfg->iterations - 1 ; i++ )
00254             {
00255                 med[i+1] = medianImage(med[i], -cfg->factor*stats->cleanstdev) ;
00256                 if (med[i+1] == NULL)
00257             {
00258                     cpl_msg_error ( _id," medianImage failed (2)" ) ;
00259                     return -1 ;
00260                 }
00261             }
00262 
00263         }
00264         else
00265         {
00266 
00267             med[0] = medianImage(threshIm, -cfg->factor) ;
00268             if ( med[0] == NULL )
00269             {
00270                 cpl_msg_error ( _id," medianImage failed (1)" ) ;
00271                 return -1 ;
00272             }
00273 
00274             for ( i = 0 ; i < cfg->iterations - 1 ; i++ )
00275             {
00276                 med[i+1] = medianImage(med[i], -cfg->factor) ;
00277                 if (med[i+1] == NULL)
00278                 {
00279                     cpl_msg_error (_id, " medianImage failed (2)" ) ;
00280                     return -1 ;
00281                 }
00282             }
00283         }
00284     }   
00285     else if (cfg->methodInd == 2)
00286     {
00287 
00288         med[0] = absDistImage(threshIm, -cfg->factor) ;
00289         if ( med[0] == NULL )
00290         {
00291             cpl_msg_error (_id, " absDistImage failed (1)" ) ;
00292             return -1 ;
00293         }
00294         for ( i = 0 ; i < cfg->iterations -1 ; i++ )
00295         {
00296             med[i+1] = absDistImage(med[i], -cfg->factor) ;
00297             if (med[i+1] == NULL)
00298             {
00299                 cpl_msg_error (_id, " absDistImage failed (2)" ) ;
00300                 return -1 ;
00301             }
00302         }
00303     }
00304     else if (cfg->methodInd == 3)
00305     {
00306         med[0] = meanImageInSpec(threshIm, -cfg->factor*stats->cleanstdev) ;
00307         if ( med[0] == NULL )
00308         {
00309             cpl_msg_error (_id, " meanImageInSpec failed (1)" ) ;
00310             return -1 ;
00311         }
00312         for ( i = 0 ; i < cfg->iterations - 1 ; i++ )
00313         {
00314             med[i+1] = meanImageInSpec(med[i], -cfg->factor*stats->cleanstdev) ;
00315             if (med[i+1] == NULL)
00316             {
00317                 cpl_msg_error (_id, " meanImageInSpec failed (2)" ) ;
00318                 return -1 ;
00319             } 
00320         }
00321     }
00322     else if (cfg->methodInd == 4)
00323     {
00324         half_box_size = (cfg->urx - cfg->llx) / 2 ;
00325         med[0] = localMedianImage(threshIm, 
00326                                   -cfg->factor, 
00327                                    cfg->loReject, 
00328                                    cfg->hiReject, 
00329                                    half_box_size) ;
00330 
00331         if ( med[0] == NULL )
00332         {
00333             cpl_msg_error ( _id," localMedianImage failed (1)" ) ;
00334             return -1 ;
00335         }
00336         for ( i = 0 ; i < cfg->iterations - 1 ; i++ )
00337         {
00338             med[i+1] = localMedianImage(med[i], 
00339                                         -cfg->factor, 
00340                                          cfg->loReject, 
00341                                          cfg->hiReject, 
00342                                          half_box_size) ;
00343 
00344             if (med[i+1] == NULL)
00345             {
00346                 cpl_msg_error (_id, " localMedianImage failed (2)" ) ;
00347                 return -1 ;
00348             }
00349         }
00350     }
00351     else
00352     {
00353         cpl_msg_error (_id, " wrong indicator methodInd in the .ini file!" ) ;
00354         return -1 ;
00355     }
00356     
00357     /* compare the filtered image with the input image */
00358     compImage = compareImages(threshIm, med[cfg->iterations - 1], medImage) ;
00359     if (compImage == NULL )
00360     {
00361         cpl_msg_error ( _id," compareImages failed" ) ;
00362         return -1 ;
00363     }
00364     
00365     /* generate the bad pixel mask */
00366     cpl_msg_info(_id,"Generates bad pixel map");
00367     maskImage = promoteImageToMask( compImage, &n ) ;
00368     if (maskImage == NULL)
00369     {
00370         cpl_msg_error (_id, " error in promoteImageToMask" ) ;
00371         return -1 ;
00372     }
00373     bp_img=cpl_image_wrap_float(maskImage->lx,
00374                                 maskImage->ly,
00375                                 maskImage->data);
00376 
00377   
00378 
00379 
00380     cpl_msg_info (_id, "no of bad pixels:  %d\n", n ) ;
00381 
00382     qclog_tbl = sinfoni_qclog_init(2);
00383     key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00384     p = cpl_parameterlist_find(config, "sinfoni.bp.method");
00385     sprintf(key_value, "%s", cpl_parameter_get_string(p));
00386  
00387     sinfoni_qclog_add(qclog_tbl,0,(char*) "QC BP-MAP METHOD",(char*) "CPL_TYPE_STRING",
00388               key_value,(char*) "BP search method");
00389 
00390     sprintf(key_value,"%d",n);
00391     sinfoni_qclog_add(qclog_tbl,1, (char*) "QC BP-MAP NBADPIX",(char*) "CPL_TYPE_INT",
00392               key_value, (char*) "No of bad pixels");
00393 
00394     if(-1 == sinfoni_pro_save_ima(bp_img,raw,sof,cfg->outName,
00395          procatg,qclog_tbl,_id,config)) {
00396          cpl_msg_error(_id,"cannot save ima %s", cfg->outName);
00397          cpl_free(key_value);
00398          for ( i = 0 ; i < cfg->iterations ; i++ )
00399              destroy_image(med[i]) ;
00400          destroy_image(medImage) ;
00401          destroy_image(medIm) ;
00402          destroy_image(compImage) ;
00403          destroy_image(maskImage) ;
00404          destroy_image(threshIm) ;
00405          if (cfg->threshInd == 1) destroy_image(colImage) ;
00406          if (stats) cpl_free(stats) ;
00407          if(med) cpl_free(med) ;
00408          badsky_free(cfg) ;
00409     
00410          return -1;
00411     }
00412     cpl_table_delete(qclog_tbl);
00413          cpl_free(key_value);
00414 
00415     for ( i = 0 ; i < cfg->iterations ; i++ )
00416         destroy_image(med[i]) ;
00417     destroy_image(medImage) ;
00418     destroy_image(medIm) ;
00419     destroy_image(compImage) ;
00420     destroy_image(maskImage) ;
00421     destroy_image(threshIm) ;
00422     if (cfg->threshInd == 1) destroy_image(colImage) ;
00423     if (stats) cpl_free(stats) ;
00424     if(med) cpl_free(med) ;
00425     badsky_free(cfg) ;
00426     
00427     
00428     return 0 ;
00429 }
00430 

Generated on Wed Oct 26 13:08:51 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001