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 ======================================= */
00044 /* =========================================================== */
00045 
00046 void spatial_filter_med_cube(TIGERfile *incube, TIGERfile *medcube,
00047                              int ring, int iter)
00048 {
00049   SPECTRUM *signal, *noise, medspec, medvar;
00050   char text[80];
00051   int *no, *idx, *nghb, *work;
00052   int i, j, k, l, nmed, nnghb, icen, npts;
00053   float *xl, *yl, delta[3];
00054   float d2max, d2, lbda;
00055   double *norm, *val, *var;
00056   double dtmp;
00057 
00058   xl  = (float *)malloc(incube->nbspec*sizeof(float)); /* Input coordinates */
00059   yl  = (float *)malloc(incube->nbspec*sizeof(float));
00060   no  =   (int *)malloc(incube->nbspec*sizeof(int));
00061   nghb=   (int *)malloc(incube->nbspec*sizeof(int));   /* Neighbor index */
00062   idx =   (int *)malloc(incube->nbspec*sizeof(int));   /* Table row (unused) */
00063 
00064   /* Get cube spatial coordinates */
00065 
00066   get_lenses_coord(incube,"X0","Y0",no,xl,yl,idx);
00067 
00068   if (RD_desc(incube,LENSIZE3,FLOAT,3,delta) < 0) {
00069     print_error("Cannot find descriptor %s in %s",LENSIZE3,incube->name);
00070 
00071     print_warning("DEPRECATED:\n"
00072                   "Previous versions of extract_spec (bef. 10/11/04,\n v1.31)\n "
00073                   "were writing an invalidly-long keyword LENSSIZEn (n=1...3).\n"
00074                   "See if you can find it back and rename it to LENSIZEn.");
00075     exit_session(ERR_BAD_DESC);
00076   }
00077   d2max = SQ((ring + 0.25)*delta[0]); /* WHY 0.25??? */
00078 
00079   /* Main loop: look for neighbors and median-filter */
00080 
00081   sprintf(text, "Iter %d - Spatial median filtering", iter);
00082   reset_print_progress();
00083   for (i=0; i<incube->nbspec; i++) {
00084     if (!DEBUG) print_progress(text,(i+1.)/incube->nbspec*100.,1.);
00085 
00086     /* Neighbors of current lens */
00087     for (nnghb=0, j=0; j<incube->nbspec; j++) {
00088       d2 = SQ(xl[j] - xl[i]) + SQ(yl[j] - yl[i]);
00089       if (d2 < d2max) {
00090         if (j == i) icen = nnghb;            /* Current lens index */
00091         nghb[nnghb++] = no[j]; 
00092       }
00093     }
00094 
00095     /* Extract neighbor spectra and compute norm */
00096     signal = (SPECTRUM *)malloc(nnghb*sizeof(SPECTRUM));
00097     noise  = (SPECTRUM *)malloc(nnghb*sizeof(SPECTRUM));
00098     norm   =   (double *)malloc(nnghb*sizeof(double));
00099 
00100     for (j=0; j<nnghb; j++) { 
00101       get_tiger_spec(incube,&(signal[j]),&(noise[j]),nghb[j]);
00102       /* Norm: use spectral mean (could be median) */
00103       // for (norm[j]=0, l=0; l<signal[j].npts; l++)
00104       //   norm[j] += RD_spec(&(signal[j]),l);
00105       // norm[j] /= signal[j].npts;
00106       /* DIRTY: Use positive-values only... */
00107       for (norm[j]=0, npts=0, l=0; l<signal[j].npts; l++)
00108         if (RD_spec(&(signal[j]),l) > 0) {
00109           norm[j] += RD_spec(&(signal[j]),l);
00110           npts ++;
00111         }
00112       if (npts) norm[j] /= npts;
00113       else print_warning("No positive pixels in spectrum #%d",nghb[j]);
00114       /* BEWARE: Noise spectrum may have previously been cosmic-flaggued */
00115       interpolate_noise(&(noise[j]));
00116     }
00117 
00118     val  = (double *)malloc(nnghb*sizeof(double));
00119     var  = (double *)malloc(nnghb*sizeof(double));
00120     work =    (int *)malloc(nnghb*sizeof(int));
00121 
00122     /* Creation of the median spectrum and variance spectrum */
00123     init_new_tiger_spec(medcube, &medspec,signal[icen].npts,signal[icen].start);
00124     init_new_tiger_spec(medcube, &medvar, signal[icen].npts,signal[icen].start);
00125 
00126     for (j=0; j<signal[icen].npts; j++) { /* Loop on the central spectrum pixels */
00127 
00128       lbda = coord_spec(&(signal[icen]),j);
00129       
00130       /* Compute median value at given lambda of normalized neighbors */
00131       for (nmed=0, k=0; k<nnghb; k++) {
00132         if (in_spectrum(&(signal[k]),lbda)) {
00133           l = pixel_spec(&(signal[k]),lbda);
00134           dtmp = norm[icen]/norm[k];
00135           val[nmed] = RD_spec(&(signal[k]),l)*dtmp;
00136           var[nmed] = RD_spec(&(noise[k]), l)*SQ(dtmp);
00137           nmed++;
00138         }
00139       }
00140 
00141       if (nmed>=3) {                        /* Compute median */
00142         WR_spec(&medspec,j,median(val,nmed,work));
00143         WR_spec(&medvar, j,median(var,nmed,work));
00144         /* Pour la variance on calcule la mediane des spectres de
00145            variance. Ce n'est probablement pas correct mais ca suffit
00146            pour ce qu'on veut en faire...(EP) */
00147       }
00148       else {                                /* Store central values */
00149         WR_spec(&medspec,j,RD_spec(&(signal[icen]),j) );
00150         WR_spec(&medvar, j,RD_spec(&(noise[icen]),j)  );
00151       }
00152 
00153     } /* Loop on the central spectrum pixels */
00154 
00155     put_tiger_spec(medcube,&medspec,&medvar,no[i]);
00156 
00157     free(val); free(var); free(work); free(norm);
00158     for (j=0; j<nnghb; j++) { 
00159       free_spec_mem(&(signal[j]));
00160       free_spec_mem(&(noise[j]));
00161     }
00162     free(signal); free(noise); 
00163 
00164   } /* Spatial loop */
00165 
00166   free(xl); free(yl); free(no); free(idx); free(nghb);
00167 }
00168 
00169 /* === Doxygen Comment ======================================= */
00183 /* =========================================================== */
00184 
00185 void spectro_spatial_filter(TIGERfile *incube, TIGERfile *medcube,
00186                             TIGERfile *rescube, int radius, int iter)
00187 {
00188   SPECTRUM medspec,inspec,resspec;
00189   char text[80];
00190   int *no;
00191   int i,j;
00192   float *diffmed, *f_diffmed;
00193 
00194   no = (int *) malloc(medcube->nbspec * sizeof(int));
00195   get_lenses_no(medcube,no);
00196 
00197   sprintf(text, "Iter %d - Spectro-spatial filtering", iter);
00198   reset_print_progress();
00199   for (i=0; i<medcube->nbspec; i++) {
00200     print_progress(text,(i+1.)/medcube->nbspec*100.,1.);
00201 
00202     get_tiger_spec(incube, &inspec, NULL,no[i]);  /* Input spectrum */
00203     get_tiger_spec(medcube,&medspec,NULL,no[i]);  /* Spat.med. spectrum */
00204 
00205     init_new_tiger_spec(rescube,&resspec,inspec.npts,inspec.start); /* Residues */
00206 
00207     diffmed   = (float*)malloc(inspec.npts*sizeof(float));
00208     f_diffmed = (float*)malloc(inspec.npts*sizeof(float));
00209 
00210     /* diffmed = spectrum - spatial median */
00211     for (j=0; j<medspec.npts; j++) 
00212       diffmed[j] = RD_spec(&inspec,j) - RD_spec(&medspec,j);
00213 
00214     /* Spectral median filter of diffmed */
00215     filter_med_array(diffmed,inspec.npts,radius,f_diffmed);
00216 
00217     /* Residual = diffmed - f_diffmed */
00218     for (j=0; j<resspec.npts; j++)
00219       WR_spec(&resspec,j, diffmed[j] - f_diffmed[j]);
00220 
00221     put_tiger_spec(rescube,&resspec,NULL,no[i]);
00222 
00223     free(diffmed); free(f_diffmed);
00224     free_spec_mem(&medspec); free_spec_mem(&inspec);
00225   }
00226 
00227   free(no);
00228 }
00229 
00230 /* === Doxygen Comment ======================================= */
00249 /* =========================================================== */
00250 
00251 int Nuke_cosmic(SPECTRUM *signal, SPECTRUM *noise,
00252                 SPECTRUM *medspec, SPECTRUM *medvar,
00253                 SPECTRUM *resspec, float sigclip, float threshold, int pixels)
00254 {
00255   int *cosmic;
00256   int i,j,ncosm;
00257   double norm, mean_signal, mean_median;
00258 
00259   cosmic = (int *)calloc(signal->npts,sizeof(int));
00260 
00261   /* 1st-pass: detect cosmics (sigma-clipping + thresholding) */
00262   for (ncosm=0, i=0; i<signal->npts; i++) {
00263     if ( RD_spec(resspec,i) > sigclip*sqrt(RD_spec(medvar,i)) && 
00264          (RD_spec(medspec,i) < 0 || 
00265           RD_spec(resspec,i) > threshold*RD_spec(medspec,i)) ) {
00266       cosmic[i] = 1;                     /* There's a cosmic there! */
00267       ncosm++;
00268     }
00269   }
00270 
00271   /* 2nd pass: flag neighbor pixels affected by cosmics */
00272   for (i=0; i<signal->npts; i++) {
00273     if (cosmic[i] == 1) {
00274       for (j=MAX(0,i-pixels); j<=MIN(signal->npts-1,i+pixels); j++) 
00275         if (!cosmic[j]) cosmic[j] = 2;   /* Flag previously clean px only */
00276     }
00277   }
00278 
00279   /* Compute normalization factor between signal and spatial median */
00280   mean_median = 0;
00281   mean_signal = 0;
00282   for (i=0; i<signal->npts; i++) { 
00283     if (!cosmic[i]) {           /* Discard cosmics from mean */
00284       mean_signal += RD_spec(signal,i);
00285       mean_median += RD_spec(medspec,i);
00286     }
00287   }
00288   norm = mean_signal / mean_median;
00289 
00290   /* Nuke cosmics and neighbor pixels */
00291   for (i=0; i<signal->npts; i++) {
00292     if (cosmic[i]) {         /* Replace signal w/ normalized spatial median */
00293       WR_spec(signal,i,RD_spec(medspec,i)*norm);
00294       if (noise != NULL) WR_spec(noise,i,HITVAR); /* Flag noise spectrum */
00295     }
00296     else WR_spec(signal,i,RD_spec(signal,i));     /* No cosmic: just copy */
00297   }
00298 
00299   free(cosmic); 
00300   
00301   return(ncosm);
00302 }
00303 
00304 /* === Doxygen Comment ======================================= */
00317 /* =========================================================== */
00318 
00319 int Remove_cosmic(TIGERfile *incube, TIGERfile *medcube, TIGERfile *rescube,
00320                   float sigclip, float threshold, int pixels, int iter)
00321 {
00322   SPECTRUM signal, medspec, medvar, resspec, noise;
00323   char text[80];
00324   int *no;
00325   int i,n,ncosmic;
00326 
00327   no = (int *)malloc(incube->nbspec * sizeof(int));
00328   get_lenses_no(incube,no);
00329 
00330   sprintf(text, "Iter %d - Removing cosmics", iter);
00331   reset_print_progress();
00332   for (ncosmic=0, i=0; i<incube->nbspec; i++) {
00333     print_progress(text,(i+1.)/incube->nbspec*100.,1.);
00334 
00335     get_tiger_spec(incube, &signal, &noise, no[i]); /* Input */
00336     get_tiger_spec(medcube,&medspec,&medvar,no[i]); /* Median */
00337     get_tiger_spec(rescube,&resspec,NULL,   no[i]); /* Residual */
00338 
00339     /* Cosmic cleaning of the spectrum */  
00340     n = Nuke_cosmic(&signal, &noise, &medspec,&medvar, &resspec,
00341                     sigclip,threshold,pixels); 
00342     ncosmic += n;                      /* Total nb of cosmics in input cube */
00343 
00344     put_tiger_spec(incube,&signal,&noise,no[i]);
00345 
00346     free_spec_mem(&medspec); free_spec_mem(&medvar); free_spec_mem(&resspec);
00347   }
00348 
00349   free(no);
00350 
00351   return(ncosmic);
00352 }
00353 
00354 /* === Doxygen Comment ======================================= */
00402 /* =========================================================== */
00403 
00404 int main(int argc, char **argv)
00405 {
00406   TIGERfile incube, outcube, medcube, rescube, dbgmedcube, dbgrescube;
00407   SPECTRUM spec, noise;
00408 
00409   char **label, **param;
00410   char *inname, *outname, *medname, *resname, *path, *name;
00411   char dbgmedname[lg_name+1], dbgresname[lg_name+1];
00412 
00413   int *no;
00414   int i, iter, itermax, ring, pixels, ncosmic, ncosmic_tot,rmed; 
00415   int fclass;
00416 
00417   float sigclip,threshold;
00418 
00419   set_purpose("Detection & suppression of cosmics in datacube");
00420   set_arglist("-in none -out none -sigclip|N 10 -ring 1 -itermax 5 -rspec 7 "
00421               "-pixels 2 -threshold 1.5 "
00422               "-dbgmed dbg_cosm_med -dbgres dbg_cosm_res");
00423 
00424   init_snifs("$Name:  $","$Revision: 1.11 $");
00425   init_session(argv, argc, &label, &param);
00426   
00427   if (DEBUG) {
00428     print_msg("$Id: remcosmic.c,v 1.11 2005/09/15 21:43:14 ycopin Exp $");
00429     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00430   }
00431 
00432   inname  = param[0];            /* Input datacube */
00433   outname = param[1];            /* Output datacube */
00434   get_argval(2,"%f",&sigclip);   /* Detection threshold sigma-clipping */
00435   get_argval(3,"%d",&ring);      /* Nb of rings for neighbor search */
00436   get_argval(4,"%d",&itermax);   /* Max nb of iterations */
00437   get_argval(5,"%d",&rmed);      /* Spectral median filtering radius */
00438   get_argval(6,"%d",&pixels);    /* Spoiled pixels on each side of cosmic */
00439   get_argval(7,"%f",&threshold); /* Relative threshold of detection */
00440 
00441   if (DEBUG) {
00442     medname = param[8];          /* Spatial median-filtered datacube */
00443     resname = param[9];          /* Residual datacube */
00444     remove_file_extension(medname);
00445     remove_file_extension(resname);
00446   }
00447   else {
00448     path = get_path(outname);                /* Handle path */
00449     name = strdup(inname);
00450     remove_path(name);                       /* Keep basic name only */
00451     remove_file_extension(name);
00452     medname = (char *)malloc((lg_name+1)*sizeof(char));
00453     resname = (char *)malloc((lg_name+1)*sizeof(char));
00454     sprintf(medname,"%s/dbg_cosm_med_%s",path,name);
00455     sprintf(resname,"%s/dbg_cosm_res_%s",path,name);
00456     free(name);
00457   }
00458 
00459   print_msg("o Input datacube: '%s'",inname);
00460   print_msg("o Output datacube: '%s'",outname);
00461   print_msg("o Sigma-clipping detection threshold: %.1f",sigclip);
00462   print_msg("o Relative threshold: %4.1f",threshold);
00463   print_msg("o Maximum nb of iterations: %d",itermax);
00464   print_msg("o Cosmic spoiled radius: %d px",pixels);
00465   print_msg("o Spectral median-filtering radius: %d px",rmed);
00466   print_msg("o Spatial median-filtering radius: %d lenses",ring);
00467   if (DEBUG) {
00468     print_msg("o Spatial median-filtered datacube: '%s'",medname);
00469     print_msg("o Residual datacube: '%s'",resname);
00470   }
00471 
00472   /* Open input datacube and test FCLASS */
00473 
00474   print_msg("Opening input datacube %s", inname);
00475   if (open_tiger_frame(&incube, inname, "I") < 0){
00476     print_error("Cannot open datacube %s",inname);
00477     exit_session(ERR_OPEN);
00478   }
00479   if (!exist_noise(&incube)) {
00480     print_error("Input datacube does not have a variance extension.");
00481     close_tiger_frame(&incube);
00482     exit_session(ERR_BAD_TYPE);    
00483   }
00484 
00485   switch ((fclass = read_file_class(&incube))) {
00486   case SUPER_CLASS:
00487   case FLA_OBJ_CUBE:
00488   case FLA_CON_CUBE:
00489   case FLA_DOM_CUBE:
00490   case FLA_SKY_CUBE: break;
00491   default :
00492     print_error("Wrong fclass (%d) for %s (should be "
00493                 "LFFF Datacube)",fclass,incube.name);
00494     print_warning("DEBUG: Pb with datacube FCLASS, keep going...");
00495     //      exit_session(ERR_BAD_TYPE);
00496     break;
00497   }
00498 
00499   /* Create output datacube */
00500 
00501   print_msg("Creating output datacube %s", outname);  
00502   if (create_tiger_frame(&outcube,outname, 1,-1, incube.step,incube.data_type, 
00503                          incube.table_name, incube.ident, incube.cunit) < 0) {
00504     close_tiger_frame(&incube);
00505     print_error("Cannot create datacube %s",outname);
00506     exit_session(ERR_CREAT);
00507   }
00508 
00509   /* Create intermediate datacubes */
00510 
00511   disable_erase_flag();
00512   if (create_tiger_frame(&medcube,medname,-1,-1, incube.step,incube.data_type, 
00513                          incube.table_name, incube.ident, incube.cunit) < 0 ||
00514       create_tiger_frame(&rescube,resname,-1,-1, incube.step,incube.data_type, 
00515                          incube.table_name, incube.ident, incube.cunit) < 0) {
00516     close_tiger_frame(&incube);
00517     close_tiger_frame(&outcube);
00518     print_error("Cannot create intermediate datacubes %s and/or %s",
00519                 medname,resname);
00520     exit_session(ERR_CREAT);
00521   }
00522   restore_erase_flag();
00523 
00524   /* Duplicate input datacube into output datacube */
00525 
00526   no = (int *)malloc(incube.nbspec*sizeof(int));
00527   get_lenses_no(&incube,no);
00528 
00529   for (i=0; i<incube.nbspec; i++) {
00530     get_tiger_spec(&incube,  &spec, &noise, no[i]);
00531     put_tiger_spec(&outcube, &spec, &noise, no[i]);
00532   }
00533 
00534   CP_non_std_desc(&incube,&outcube);
00535   switch (fclass) {
00536   case FLA_OBJ_CUBE: write_file_class(&outcube, COS_OBJ_CUBE); break;
00537   case FLA_CON_CUBE: write_file_class(&outcube, COS_CON_CUBE); break;
00538   case FLA_DOM_CUBE: write_file_class(&outcube, COS_DOM_CUBE); break;
00539   case FLA_SKY_CUBE: write_file_class(&outcube, COS_SKY_CUBE); break;
00540   }
00541   
00542   /* Cosmic cleaning (up to itermax) */
00543 
00544   ncosmic_tot = 0;              /* Total nb of cosmics */
00545   ncosmic = 1;                  /* Nb of cosmic detected during last iter. */
00546   iter = 1;                     /* Iteration */
00547   while (ncosmic!=0 && iter<=itermax) {
00548     /* Spatial median-filtering */
00549     spatial_filter_med_cube(&outcube,&medcube,ring,iter); 
00550 
00551     /* Spectro-spatial residuals */
00552     spectro_spatial_filter(&outcube,&medcube,&rescube,rmed,iter); 
00553 
00554     /* Cosmics cleaning */
00555     ncosmic = Remove_cosmic(&outcube,&medcube,&rescube,sigclip,threshold,
00556                             pixels,iter);
00557 
00558     ncosmic_tot += ncosmic;
00559     print_msg("Iter %d - %d cosmics found", iter, ncosmic);
00560 
00561     if (DEBUG) {
00562       /* Store intermediate datacubes */
00563       disable_erase_flag();
00564       sprintf(dbgmedname,"%s_%d",medname,iter);
00565       sprintf(dbgresname,"%s_%d",resname,iter);
00566       print_msg("DEBUG: Saving intermediate datacubes %s and %s (iter %d)",
00567                 dbgmedname, dbgresname, iter);
00568       if (create_tiger_frame(&dbgmedcube,dbgmedname,-1,-1, 
00569                              incube.step,incube.data_type,incube.table_name,
00570                              incube.ident,incube.cunit) < 0 ||
00571           create_tiger_frame(&dbgrescube,dbgresname,-1,-1, 
00572                              incube.step,incube.data_type,incube.table_name,
00573                              incube.ident,incube.cunit) < 0) {
00574         print_error("Cannot create debug intermediate datacubes %s and/or %s",
00575                     dbgmedname,dbgresname);
00576         exit_session(ERR_CREAT);
00577       }
00578       for (i=0; i<medcube.nbspec; i++) {
00579         get_tiger_spec(&medcube,    &spec, &noise, no[i]);
00580         put_tiger_spec(&dbgmedcube, &spec, &noise, no[i]);
00581         get_tiger_spec(&rescube,    &spec, NULL,   no[i]);
00582         put_tiger_spec(&dbgrescube, &spec, NULL,   no[i]);
00583       }
00584       write_file_class(&dbgmedcube, DEBUG_FILE); 
00585       write_file_class(&dbgrescube, DEBUG_FILE); 
00586       close_tiger_frame(&dbgmedcube);
00587       close_tiger_frame(&dbgrescube);
00588 
00589       restore_erase_flag();
00590     }
00591 
00592     iter ++;
00593   }
00594 
00595   print_msg(" -> %d cosmics found in total", ncosmic_tot);
00596   WR_desc(&outcube, COSMNB, INT, 1, &ncosmic_tot);
00597 
00598   /* Closing */
00599 
00600   close_tiger_frame(&incube);
00601   close_tiger_frame(&outcube);
00602 
00603   delete_tiger_frame(&medcube); /* They were already saved in DEBUG-mode */
00604   delete_tiger_frame(&rescube);
00605 
00606   free(no);
00607   if (!DEBUG) {
00608     free(medname);
00609     free(resname);
00610   }
00611 
00612   exit_session(OK);
00613   return(OK);
00614 }

Generated on Wed Oct 26 23:59:39 2005 for Snifs by doxygen 1.3.5