00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "bp_sky.h"
00013 #include "badsky_ini_by_cpl.h"
00014 #include "sinfoni_pro_save.h"
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
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
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
00094
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
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
00130
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
00149 destroy_cube(cubeflat) ;
00150 cpl_free(imagelist) ;
00151
00152
00153
00154
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
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
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
00239
00240
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
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
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