bp_norm.c

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

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