Main Page | Modules | Data Structures | File List | Data Fields | Globals | Related Pages

calib/source/remcosmic.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00014 /* =========================================================== */
00015 
00016 #include <IFU_io.h>
00017 #include <IFU_math.h>
00018 #include <snifs.h>
00019 #include <calib.h>
00020 
00022 static int const HITVAR=-1;
00023 
00024 /* === Doxygen Comment ======================================= */
00037 /* =========================================================== */
00038 
00039 void spatial_filter_med_cube(TIGERfile *incube, TIGERfile *medcube,
00040                              int ring, int iter)
00041 {
00042   SPECTRUM *signal, *noise, medspec, medvar;
00043   char text[80];
00044   int *no, *idx, *nghb, *work;
00045   int i, j, k, l, n, nnghb, icen;
00046   float *xl, *yl, delta[3];
00047   float d2max, d2, lbda;
00048   double *norm, *val, *var;
00049 
00050   xl  = (float *)malloc(incube->nbspec*sizeof(float)); /* Input coordinates */
00051   yl  = (float *)malloc(incube->nbspec*sizeof(float));
00052   no  =   (int *)malloc(incube->nbspec*sizeof(int));
00053   idx =   (int *)malloc(incube->nbspec*sizeof(int));
00054   nghb=   (int *)malloc(incube->nbspec*sizeof(int));   /* Neighbors */
00055 
00056   /* Get cube spatial coordinates */
00057 
00058   get_lenses_coord(incube,"X0","Y0",no,xl,yl,idx);
00059 
00060   if (RD_desc(incube,LENSIZE3,FLOAT,3,delta) < 0) {
00061     print_error("Cannot find descriptor %s in intput datacube",LENSIZE3);
00062 
00063     print_warning("Previous versions of extract_spec (bef. 10/11/04, v.1.31)\n"
00064                   "were writing an invalidly-long keyword LENSSIZEn (n=1...3).\n"
00065                   "See if you can find it back and rename it to LENSIZEn.\n");
00066 
00067     exit_session(ERR_BAD_DESC);
00068   }
00069   d2max = SQ((ring + 0.25)*delta[0]); /* WHY 0.25??? */
00070 
00071   /* Main loop: look for neighbors and median-filter */
00072 
00073   sprintf(text, "Iter %d - Spatial median filtering", iter);
00074   reset_print_progress();
00075   for (i=0; i<incube->nbspec; i++) {
00076     if (!DEBUG) print_progress(text,(i+1.)/incube->nbspec*100.,1.0);
00077 
00078     /* Neighbors of current lens */
00079     for (nnghb=0, j=0; j<incube->nbspec; j++) {
00080       d2 = SQ(xl[j] - xl[i]) + SQ(yl[j] - yl[i]);
00081       if (d2 < d2max) {
00082         if (j == i) icen = nnghb; /* Current lens index */
00083         nghb[nnghb++] = no[j]; 
00084       }
00085     }
00086 
00087     /* Extract neighbor spectra and compute norm */
00088     signal = (SPECTRUM *)malloc(nnghb*sizeof(SPECTRUM));
00089     noise  = (SPECTRUM *)malloc(nnghb*sizeof(SPECTRUM));
00090     norm   =   (double *)malloc(nnghb*sizeof(double));
00091 
00092     for (j=0; j<nnghb; j++) { 
00093       get_tiger_spec(incube,&(signal[j]),&(noise[j]),nghb[j]);
00094       for (norm[j]=0, l=0; l<signal[j].npts; l++)
00095         norm[j] += RD_spec(&(signal[j]),l);
00096       norm[j] /= signal[j].npts;
00097     }
00098 
00099     val  = (double *)malloc(nnghb*sizeof(double));
00100     var  = (double *)malloc(nnghb*sizeof(double));
00101     work =    (int *)malloc(nnghb*sizeof(int));
00102 
00103     /* Creation of the median spectrum and variance spectrum */
00104     init_new_tiger_spec(medcube, &medspec,signal[icen].npts,signal[icen].start);
00105     init_new_tiger_spec(medcube, &medvar, signal[icen].npts,signal[icen].start);
00106 
00107     for (j=0; j<signal[icen].npts; j++) { /* Loop on the central spectrum pixels */
00108 
00109       lbda = coord_spec(&(signal[icen]),j);
00110       
00111       /* Compute median value at given lambda of normalized neighbors */
00112       for (n=0, k=0; k<nnghb; k++) {
00113         if (in_spectrum(&(signal[k]),lbda)) {
00114           l = pixel_spec(&(signal[k]),lbda);
00115           val[n] = RD_spec(&(signal[k]),l)/norm[k] * norm[icen];
00116           var[n] = RD_spec(&(noise[k]), l)/norm[k] * norm[icen];
00117           n++;
00118         }
00119       }
00120 
00121       if (n >= 3) {                          /* Compute median */
00122         WR_spec(&medspec,j,median(val,n,work));
00123         WR_spec(&medvar, j,median(var,n,work));
00124         /* Pour la variance on calcule la mediane des spectres de
00125            variance. Ce n'est probablement pas correct mais ca suffit
00126            pour ce qu'on veut en faire...(EP) */
00127       }
00128       else {                                 /* Store central values */
00129         WR_spec(&medspec,j, RD_spec(&(signal[icen]),j) );
00130         WR_spec(&medvar, j, RD_spec(&(noise[icen]),j)  );
00131       }
00132 
00133     } /* Loop on the central spectrum pixels */
00134 
00135     put_tiger_spec(medcube,&medspec,&medvar,no[i]);
00136 
00137     free(val); free(var); free(work); 
00138     free(signal); free(noise); free(norm);
00139 
00140   } /* Spatial loop */
00141 
00142   free(xl); free(yl); free(no); free(idx); free(nghb);
00143 }
00144 
00145 /* === Doxygen Comment ======================================= */
00159 /* =========================================================== */
00160 
00161 void spectro_spatial_filter(TIGERfile *incube, TIGERfile *medcube,
00162                             TIGERfile *rescube, int radius, int iter)
00163 {
00164   SPECTRUM sm,sc,sr;
00165   char text[80];
00166   int *no;
00167   int i,j;
00168   float *diffmed, *f_diffmed;
00169 
00170   no = (int *) malloc(medcube->nbspec * sizeof(int));
00171   get_lenses_no(medcube,no);
00172 
00173   sprintf(text, "Iter %d - Spectro-spatial filtering", iter);
00174   reset_print_progress();
00175   for (i=0; i<medcube->nbspec; i++) {
00176     print_progress(text,(i+1.)/medcube->nbspec*100.,1.);
00177 
00178     get_tiger_spec(incube, &sc,NULL,no[i]);  /* Input spectrum */
00179     get_tiger_spec(medcube,&sm,NULL,no[i]);  /* Spat.med. spectrum */
00180 
00181     init_new_tiger_spec(rescube,&sr,sc.npts,sc.start); /* Residual spec. */
00182 
00183     diffmed   = (float*)malloc(sc.npts*sizeof(float));
00184     f_diffmed = (float*)malloc(sc.npts*sizeof(float));
00185 
00186     /* diffmed = spectrum - spatial median */
00187     for (j=0; j<sm.npts; j++) 
00188       diffmed[j] = RD_spec(&sc,j) - RD_spec(&sm,j);
00189 
00190     /* Spectral median filter of diffmed */
00191     filter_med_array(diffmed,sc.npts,radius,f_diffmed);
00192 
00193     /* Residual = diffmed - f_diffmed */
00194     for (j=0; j<sr.npts; j++)
00195       WR_spec(&sr,j, diffmed[j] - f_diffmed[j]);
00196 
00197     put_tiger_spec(rescube,&sr,NULL,no[i]);
00198 
00199     free(diffmed); free(f_diffmed);
00200     free_spec_mem(&sm); free_spec_mem(&sc);
00201   }
00202 
00203   free(no);
00204 }
00205 
00206 /* === Doxygen Comment ======================================= */
00223 /* =========================================================== */
00224 
00225 int Nuke_cosmic(SPECTRUM *signal, SPECTRUM *noise,
00226                 SPECTRUM *medspec, SPECTRUM *medvar,
00227                 SPECTRUM *resspec, float sigclip, float threshold, int pixels)
00228 {
00229   int *cosmic;
00230   int i,j,ncosm;
00231   double sigma, norm, mean_signal, mean_median;
00232 
00233   cosmic = (int *)calloc(signal->npts,sizeof(int));
00234 
00235   /* 1st-pass: detect cosmics (sigma-clipping + thresholding) */
00236   for (ncosm=0, i=0; i<signal->npts; i++) {
00237     sigma = sqrt(RD_spec(medvar,i));         /* Spatial median variance */
00238     if ((RD_spec(resspec,i) > sigclip*sigma) && 
00239         (RD_spec(resspec,i)/RD_spec(medspec,i) > threshold)) {
00240       cosmic[i] = 1;                         /* There's a cosmic there! */
00241       ncosm++;
00242     }
00243   }
00244 
00245   /* 2nd pass: flag neighbor pixels affected by cosmics */
00246   for (i=0; i<signal->npts; i++) {
00247     if (cosmic[i] == 1) {
00248       for (j=MAX(0,i-pixels); j<=MIN(signal->npts-1,i+pixels); j++) 
00249         if (!cosmic[j]) cosmic[j] = 2; /* Flag previously clean px only */
00250     }
00251   }
00252 
00253   /* Compute normalization factor between signal and spatial median */
00254   mean_median = 0;
00255   mean_signal = 0;
00256   for (i=0; i<signal->npts; i++) { 
00257     if (!cosmic[i]) {           /* Discard cosmics from mean */
00258       mean_signal += RD_spec(signal,i);
00259       mean_median += RD_spec(medspec,i);
00260     }
00261   }
00262   norm = mean_signal / mean_median;
00263 
00264   /* Nuke cosmics and neighbor pixels */
00265   for (i=0; i<signal->npts; i++) {
00266     if (cosmic[i]) {         /* Replace signal w/ normalized spatial median */
00267       WR_spec(signal,i,RD_spec(medspec,i)*norm);
00268       if (noise != NULL) WR_spec(noise,i,HITVAR); /* Flag noise spectrum */
00269     }
00270     else WR_spec(signal,i,RD_spec(signal,i)); /* No cosmic: just copy */
00271   }
00272 
00273   free(cosmic); 
00274   
00275   return(ncosm);
00276 }
00277 
00278 /* === Doxygen Comment ======================================= */
00291 /* =========================================================== */
00292 
00293 int Remove_cosmic(TIGERfile *incube, TIGERfile *medcube, TIGERfile *rescube,
00294                   float sigclip, float threshold, int pixels, int iter)
00295 {
00296   SPECTRUM sc, sm, nm, sr, noise;
00297   char text[80];
00298   int *no;
00299   int i,n,ncosmic, NOISE;
00300 
00301   no = (int *)malloc(incube->nbspec * sizeof(int));
00302   get_lenses_no(incube,no);
00303   NOISE = exist_noise(incube);
00304 
00305   sprintf(text, "Iter %d - Removing cosmics", iter);
00306   reset_print_progress();
00307   for (ncosmic=0, i=0; i<incube->nbspec; i++) {
00308     print_progress(text,(i+1.)/incube->nbspec*100.,1.);
00309 
00310     get_tiger_spec(incube, &sc,(NOISE? &noise:NULL),no[i]); /* Input */
00311     get_tiger_spec(medcube,&sm,        &nm,         no[i]); /* Median */
00312     get_tiger_spec(rescube,&sr,        NULL,        no[i]); /* Residual */
00313 
00314     /* Cosmic cleaning of the spectrum */  
00315     n = Nuke_cosmic(&sc, (NOISE? &noise:NULL), &sm,&nm, &sr,
00316                     sigclip,threshold,pixels); 
00317     ncosmic += n;                      /* Total nb of cosmics in input cube */
00318 
00319     put_tiger_spec(incube,&sc,(NOISE? &noise:NULL),no[i]);
00320 
00321     free_spec_mem(&sm); free_spec_mem(&sr);
00322   }
00323 
00324   free(no);
00325 
00326   return(ncosmic);
00327 }
00328 
00329 /* === Doxygen Comment ======================================= */
00348 /* =========================================================== */
00349 
00350 int main(int argc, char **argv)
00351 {
00352   TIGERfile incube, outcube, medcube, rescube;
00353   SPECTRUM spec, noise;
00354 
00355   char **label, **param;
00356   char *inname, *outname, *medname, *resname;
00357 
00358   int *no;
00359   int i, iter, itermax, ring, pixels, ncosmic, ncosmic_tot,rmed; 
00360   int fclass, NOISE;
00361 
00362   float sigclip,threshold;
00363 
00364   set_purpose("Detection & suppression of cosmics in datacube");
00365   set_arglist("-in none -out none -sigclip|N 10 -ring 1 -itermax 5 -rspec 3 "
00366               "-pixels 2 -threshold 1.5 -dbgmed med -dbgres res");
00367 
00368   init_snifs("$Name:  $");
00369   init_session(argv, argc, &label, &param);
00370   
00371   if (DEBUG) {
00372     print_msg("$Id: remcosmic.c,v 1.7 2004/11/10 17:45:17 ycopin Exp $");
00373     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00374   }
00375 
00376   inname  = param[0];            /* Input datacube */
00377   outname = param[1];            /* Output datacube */
00378   get_argval(2,"%f",&sigclip);   /* Detection threshold sigma-clipping */
00379   get_argval(3,"%d",&ring);      /* Nb of rings for neighbor search */
00380   get_argval(4,"%d",&itermax);   /* Max nb of iterations */
00381   get_argval(5,"%d",&rmed);      /* Spectral median filtering radius */
00382   get_argval(6,"%d",&pixels);    /* Spoiled pixels on each side of cosmic */
00383   get_argval(7,"%f",&threshold); /* Relative threshold of detection */
00384 
00385   if (DEBUG) {
00386     medname = param[8];          /* Spatial median-filtered datacube */
00387     resname = param[9];          /* Residual datacube */
00388   }
00389   else {
00390     medname = (char *)malloc((lg_name+1)*sizeof(char));
00391     resname = (char *)malloc((lg_name+1)*sizeof(char));
00392     sprintf(medname,"dbgmed_%s",param[0]);
00393     remove_file_extension(medname);
00394     sprintf(resname,"dbgres_%s",param[0]);
00395     remove_file_extension(resname);
00396   }
00397 
00398   print_msg("o Input datacube: '%s'",inname);
00399   print_msg("o Output datacube: '%s'",outname);
00400   print_msg("o Sigma-clipping detection threshold: %.1f",sigclip);
00401   print_msg("o Relative threshold: %4.1f",threshold);
00402   print_msg("o Maximum nb of iterations: %d",itermax);
00403   print_msg("o Cosmic spoiled radius: %d px",pixels);
00404   print_msg("o Median-filtering spectral radius: %d px",rmed);
00405   print_msg("o Spatial radius: %d lenses",ring);
00406   if (DEBUG) {
00407     print_msg("o Spatial median-filtered datacube: '%s'",medname);
00408     print_msg("o Residual datacube: '%s'",resname);
00409   }
00410 
00411   /* Open input datacube and test FCLASS */
00412 
00413   print_msg("Opening input datacube %s", inname);
00414   if (open_tiger_frame(&incube, inname, "I") < 0){
00415     print_error("Cannot open input datacube %s",inname);
00416     exit_session(ERR_OPEN);
00417   }
00418 
00419   fclass = read_file_class(&incube);
00420   switch (fclass) {
00421   case SUPER_CLASS:
00422   case FLA_OBJ_CUBE:
00423   case FLA_CON_CUBE:
00424   case FLA_SKY_CUBE: break;
00425   default :
00426     print_warning("Wrong file class (should be Flat-fielded Datacube)");
00427     print_warning("WE SHOULD STOP HERE but we let it go for the moment");
00428     //    exit_session(ERR_BAD_TYPE);
00429     break;
00430   }
00431 
00432   /* Create output datacube */
00433 
00434   print_msg("Creating output datacube %s", outname);  
00435   if (create_tiger_frame(&outcube,outname, 1,-1, incube.step,incube.data_type, 
00436                          incube.table_name, incube.ident, incube.cunit) < 0) {
00437     close_tiger_frame(&incube);
00438     print_error("Cannot create output datacube %s",outname);
00439     exit_session(ERR_CREAT);
00440   }
00441 
00442   /* Create intermediate and debug datacubes */
00443 
00444   disable_erase_flag();
00445   if (create_tiger_frame(&medcube,medname,-1,-1, incube.step,incube.data_type, 
00446                          incube.table_name, incube.ident, incube.cunit) < 0 ||
00447       create_tiger_frame(&rescube,resname,-1,-1, incube.step,incube.data_type, 
00448                          incube.table_name, incube.ident, incube.cunit) < 0) {
00449     close_tiger_frame(&incube);
00450     close_tiger_frame(&outcube);
00451     print_error("Cannot create intermediate datacubes %s and/or %s",
00452                 medname,resname);
00453     exit_session(ERR_CREAT);
00454   }
00455   restore_erase_flag();
00456 
00457   /* Duplicate input datacube into output datacube */
00458 
00459   no = (int *)malloc(incube.nbspec*sizeof(int));
00460   get_lenses_no(&incube,no);
00461   NOISE = exist_noise(&incube);
00462 
00463   for (i=0; i<incube.nbspec; i++) {
00464     get_tiger_spec(&incube,  &spec, (NOISE? &noise:NULL), no[i]);
00465     put_tiger_spec(&outcube, &spec, (NOISE? &noise:NULL), no[i]);
00466   }
00467 
00468   CP_non_std_desc(&incube,&outcube);
00469   switch (fclass) {
00470   case FLA_OBJ_CUBE: write_file_class(&outcube, COS_OBJ_CUBE); break;
00471   case FLA_CON_CUBE: write_file_class(&outcube, COS_CON_CUBE); break;
00472   case FLA_SKY_CUBE: write_file_class(&outcube, COS_SKY_CUBE); break;
00473   }
00474   
00475   /* Cosmic cleaning (up to itermax) */
00476 
00477   ncosmic_tot = 0;              /* Total nb of cosmics */
00478   ncosmic = 1;                  /* Nb of cosmic detected during last iter. */
00479   iter = 1;                     /* Iteration */
00480   while ((ncosmic!=0) && (iter <= itermax)) {
00481     /* Spatial median-filtering */
00482     spatial_filter_med_cube(&outcube,&medcube,ring,iter); 
00483 
00484     /* Spectro-spatial residuals */
00485     spectro_spatial_filter(&outcube,&medcube,&rescube,rmed,iter); 
00486 
00487     /* Cosmics cleaning */
00488     ncosmic = Remove_cosmic(&outcube,&medcube,&rescube,sigclip,threshold,
00489                             pixels,iter);
00490 
00491     ncosmic_tot += ncosmic;
00492     print_msg("Iter %d - %d cosmics found", iter, ncosmic);
00493     iter ++;
00494   }
00495 
00496   print_msg(" -> %d cosmics found in total", ncosmic_tot);
00497   WR_desc(&outcube, "NBCOSM", INT, 1, &ncosmic_tot);
00498 
00499   /* Closing */
00500 
00501   close_tiger_frame(&incube);
00502   close_tiger_frame(&outcube);
00503   
00504   if (DEBUG){
00505     print_msg("Saving intermediate datacubes %s and %s",
00506               medcube.name, rescube.name);
00507     write_file_class(&medcube, DEBUG_FILE); close_tiger_frame(&medcube);
00508     write_file_class(&rescube, DEBUG_FILE); close_tiger_frame(&rescube);
00509   }
00510   else {
00511     delete_tiger_frame(&medcube);
00512     delete_tiger_frame(&rescube);
00513   }
00514 
00515   exit_session(OK);
00516   return(OK);
00517 }

Generated on Tue Nov 23 18:04:18 2004 for Snifs by doxygen 1.3.3