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

calib/source/apply_lfff.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00017 /* =========================================================== */
00018 
00019 #include <IFU_io.h>
00020 #include <IFU_math.h>
00021 #include <snifs.h>
00022 #include <calib.h>
00023 
00024 /* === Doxygen Comment ======================================= */
00030 /* =========================================================== */
00031 
00032 int interWspec(SPECTRUM *spec1,SPECTRUM *spec2)                 
00033 {
00034   double l1, l2;
00035 
00036   l1 = MAX(spec1->wstart,spec2->wstart);
00037   l2 = MIN(spec1->wend,spec2->wend);
00038   if (l1 >= l2) return(ERR_BAD_PARAM);
00039   spec1->iwstart = ((l1 - spec1->start) / spec1->step) + 0.5;
00040   spec1->iwend   = ((l2 - spec1->start) / spec1->step) + 0.5;
00041   spec1->wstart  = spec1->start + spec1->iwstart*spec1->step;
00042   spec1->wend    = spec1->start + spec1->iwend*spec1->step;
00043   spec2->iwstart = ((l1 - spec2->start) / spec2->step) + 0.5;
00044   spec2->iwend   = ((l2 - spec2->start) / spec2->step) + 0.5;
00045   spec2->wstart  = spec2->start + spec2->iwstart*spec2->step;
00046   spec2->wend    = spec2->start + spec2->iwend*spec2->step;
00047   return(OK);
00048 }
00049 
00050 /* === Doxygen Comment ======================================= */
00063 /* =========================================================== */
00064 
00065 void applyLfff(TIGERfile *incube, TIGERfile *fcube, TIGERfile *outcube, int COL, 
00066                int TSKY, char *colname1, char *colname2)
00067 {
00068   SPECTRUM inspec, fspec, outspec;
00069   SPECTRUM innoise,outnoise;
00070   TABLE ftable;
00071 
00072   int *nolens, *row;
00073   int col_no, col_flag,col_lb1, col_lb2;
00074   int i, j, k, k1, k2, l, npts, nrej, ntrunc, nzero, flag, nbsel, NOISE;
00075 
00076   float var, res, valff, lb1, lb2, *A, *D;
00077 
00078   nrej = 0;
00079   ntrunc = 0;
00080 
00081   nolens = (int *)malloc(incube->nbspec*sizeof(int));
00082   A    = (float *)malloc(incube->nbspec*sizeof(float));
00083   D    = (float *)malloc(incube->nbspec*sizeof(float));
00084   row  =   (int *)malloc(incube->nbspec*sizeof(int));
00085   get_lenses_coord_select(incube,"A","D",nolens,A,D,row,&nbsel);
00086 
00087   if (COL) {
00088     print_msg("Opening Flat table %s (col. %s,%s)",
00089               fcube->table_name, colname1, colname2);
00090     open_table(&ftable, fcube->table_name, "I");
00091     col_no   = get_col_ref(&ftable, "No");
00092     col_lb1  = get_col_ref(&ftable, colname1);
00093     col_lb2  = get_col_ref(&ftable, colname2);
00094     col_flag = get_col_ref(&ftable, "FLAG");
00095   }
00096 
00097   set_control_level(NONE);
00098   handle_select_flag(&ftable, 'W', NULL);
00099   NOISE = exist_noise(incube);
00100 
00101   for (i=0; i<nbsel; i++) {                  /* Loop over flat table rows */
00102     print_progress("Flat-Fielding",100.*(i+1)/incube->nbspec,1);
00103 
00104     if (NOISE) get_tiger_spec(incube, &inspec, &innoise, nolens[i]);
00105     else       get_tiger_spec(incube, &inspec, NULL, nolens[i]);
00106 
00107     if (get_tiger_spec(fcube, &fspec, NULL, nolens[i]) != 0) {
00108       print_warning("Spectra %d discarded (no flat field spectrum)",nolens[i]);
00109       free_spec_mem(&inspec);
00110       if (NOISE) free_spec_mem(&innoise);
00111       nrej++;
00112       continue;
00113     }
00114 
00115     if (COL) {
00116       l = search_in_col(&ftable, col_no, &(nolens[i])); 
00117       if (RD_tbl(&ftable, l, col_flag, &flag) == ERR_NODATA) {
00118         print_warning("Spectra %d discarded (no common pixel)",nolens[i]);
00119         free_spec_mem(&inspec);
00120         if (NOISE) free_spec_mem(&innoise);
00121         free_spec_mem(&fspec);
00122         nrej++;
00123         continue;
00124       }
00125       /* special treatment if object or (TSKY and sky) */
00126       if ((flag == 1) || (TSKY && flag == 2)) {
00127         if ((RD_tbl(&ftable, l, col_lb1, &lb1) == ERR_NODATA) ||
00128             (RD_tbl(&ftable, l, col_lb2, &lb2) == ERR_NODATA)) 
00129           npts = 0;
00130         else {
00131           k1  = pixel_spec(&fspec, lb1); 
00132           lb1 = coord_spec(&fspec, k1);
00133           fspec.iwstart = k1; 
00134           fspec.wstart  = lb1;
00135           k2  = pixel_spec(&fspec, lb2); 
00136           lb2 = coord_spec(&fspec, k2);
00137           fspec.iwend = k2; 
00138           fspec.wend  = lb2;
00139           
00140           interWspec(&inspec, &fspec);
00141           npts = inspec.iwend - inspec.iwstart + 1; 
00142         } 
00143       } 
00144       else {                                 /* usual treatment otherwise */
00145         inter_spec(&inspec, &fspec);
00146         npts = inspec.iwend - inspec.iwstart + 1; 
00147       }
00148     } /* COL */
00149     else {                                   /* usual treatment */
00150       inter_spec(&inspec, &fspec);
00151       npts = inspec.iwend - inspec.iwstart + 1; 
00152     }
00153     
00154     if (npts <= 0) {
00155       print_warning("Spectra %d discarded (no common pixel)", nolens[i]);
00156       free_spec_mem(&inspec);
00157       if (NOISE) free_spec_mem(&innoise);
00158       free_spec_mem(&fspec);
00159       nrej++;
00160       continue;
00161     }
00162     if (npts < inspec.npts) ntrunc++;           /* truncated spectra */
00163   
00164     init_new_tiger_spec(outcube,&outspec,npts,inspec.wstart);
00165     if (NOISE) init_new_tiger_spec(outcube,&outnoise,npts,inspec.wstart);
00166 
00167     /* Flat-fielding */
00168     nzero = 0;
00169     for(k=inspec.iwstart, l=fspec.iwstart, j=0; j<npts; j++, k++, l++) {
00170       valff = RD_spec(&fspec,l);
00171       if (valff == 0) {                      /* cannot flat-field */
00172         nzero++;
00173         WR_spec(&outspec,j,0);
00174         if (NOISE) WR_spec(&outnoise,j,-1);
00175       } 
00176       else {
00177         res = RD_spec(&inspec,k)/valff;
00178         WR_spec(&outspec,j,res);
00179         if (NOISE) {
00180           var = RD_spec(&innoise,k)/(valff*valff);
00181           WR_spec(&outnoise,j,var);
00182         }
00183       }
00184     }
00185     if (nzero > 0) 
00186       print_warning("%d flagged bad points in spectrum %d",nzero, nolens[i]);
00187 
00188     /* Store the result */
00189     if (NOISE) put_tiger_spec(outcube, &outspec, &outnoise, nolens[i]);
00190     else       put_tiger_spec(outcube, &outspec, NULL, nolens[i]);
00191 
00192     free_spec_mem(&inspec);
00193     if (NOISE) free_spec_mem(&innoise);
00194     free_spec_mem(&fspec); 
00195 
00196   } /* Loop over flat table rows */
00197 
00198   if (COL) close_table(&ftable);
00199 
00200   print_msg("-> %d spectra corrected (%d truncated) and %d spectra rejected", 
00201             nbsel, ntrunc, nrej);
00202 
00203   free(nolens); free(A); free(D); free(row);
00204 }
00205 
00206 /* === Doxygen Comment ======================================= */
00219 /* =========================================================== */
00220 
00221 int main(int argc,char **argv)
00222 {
00223   TIGERfile incube, fcube, outcube;
00224   char **label, **param, **colname;
00225   char *inname, *outname, *fname;
00226   int npix = -1, class, COL, TSKY;
00227   double start=0,end=0;
00228 
00229   set_arglist("-in none -lfff none -out none -col null -tsky");
00230   set_purpose("Low-frequency spectro-spatial flat-field (application)");
00231   
00232   init_snifs("$Name:  $","$Revision: 1.11 $");
00233   init_session(argv, argc, &label, &param);
00234 
00235   if (DEBUG) {
00236     print_msg("$Id: apply_lfff.c,v 1.11 2005/10/24 10:39:50 ycopin Exp $");
00237     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00238   }
00239 
00240   /* Inputs */
00241 
00242   inname = param[0];
00243   fname  = param[1];
00244   outname= param[2];
00245 
00246   print_msg("Datacubes Input: %s Output: %s", inname, outname);
00247   print_msg("References Flatfield Datacube: %s", fname);
00248 
00249   colname    = (char **)malloc(2*sizeof(char *));
00250   colname[0] = (char *)malloc((lg_label+1)*sizeof(char));
00251   colname[1] = (char *)malloc((lg_label+1)*sizeof(char));     
00252 
00253   if ((COL = is_set(param[3]))) {
00254     if (decode_argval_char(param[3], colname) != 2) {
00255       print_error("Error in input parameter -col col(start),col(end) (got %s)", 
00256                   param[3]);
00257       exit_session(ERR_BAD_PARAM);
00258     }
00259     if ((TSKY = is_true(param[4]))) 
00260       print_msg("Object AND Sky spectra truncated according to columns %s,%s",
00261                 colname[0], colname[1]);
00262     else
00263       print_msg("Only Object spectra truncated according to columns %s,%s",
00264                 colname[0], colname[1]);
00265   }
00266   
00267   /* Open input datacube and check file class */
00268   print_msg("Opening input datacube %s",inname);
00269   if (open_tiger_frame(&incube,inname,"I") < 0){
00270     print_error("Cannot open datacube %s",inname);
00271     exit_session(ERR_OPEN);
00272   }
00273 
00274   switch ((class = read_file_class(&incube))) {
00275   case SUPER_CLASS:
00276   case WAV_OBJ_CUBE: 
00277   case WAV_SKY_CUBE:
00278   case WAV_DOM_CUBE:
00279   case WAV_CON_CUBE: break;
00280   default:
00281     print_error("Wrong fclass (%d) for %s (should be lambda-calibrated datacube)",
00282                 class,inname);
00283     exit_session(ERR_BAD_PARAM);
00284   }
00285 
00286   /* Open flat datacube and check file class */
00287   print_msg("Opening flat-field datacube %s",inname);
00288   if (open_tiger_frame(&fcube, fname, "I") < 0){
00289     close_tiger_frame(&incube);
00290     print_error("Cannot open datacube %s",fname);
00291     exit_session(ERR_OPEN);
00292   }
00293 
00294   switch (read_file_class(&fcube)) {
00295   case SUPER_CLASS:
00296   case FLAT_CUBE: break;
00297   default:
00298     print_error("Wrong fclass for %s (should be Flat-field datacube)",fname);
00299     exit_session(ERR_BAD_PARAM);
00300   }
00301 
00302   if (ABS(incube.step - fcube.step) > 0.0000001){ /* strange test... */
00303     close_tiger_frame(&incube);
00304     close_tiger_frame(&fcube);
00305     print_error("Input and Flat datacubes have different steps");
00306     exit_session(ERR_BAD_PARAM);
00307   }
00308 
00309   if (has_common_bounds(&incube)) get_common_param(&incube,&npix,&start,&end);
00310   else {
00311     npix = -1;
00312     start = -1;
00313   }
00314 
00315   if (create_tiger_frame(&outcube,outname,npix,start,incube.step,incube.data_type,
00316                          incube.table_name,incube.ident,incube.cunit) < 0){
00317     close_tiger_frame(&incube);
00318     close_tiger_frame(&fcube);
00319     print_error("Cannot create datacube %s",outname);
00320     exit_session(ERR_OPEN);
00321   }
00322 
00323   applyLfff(&incube, &fcube, &outcube, COL, TSKY, colname[0], colname[1]);
00324 
00325   /* write output datacube descriptors */
00326   CP_non_std_desc(&incube, &outcube);  
00327   WR_desc(&outcube, LFFFCUBE,  CHAR, lg_name+1, fcube.name);
00328   switch (class) {
00329   case WAV_OBJ_CUBE: write_file_class(&outcube, FLA_OBJ_CUBE); break;
00330   case WAV_SKY_CUBE: write_file_class(&outcube, FLA_SKY_CUBE); break;
00331   case WAV_DOM_CUBE: write_file_class(&outcube, FLA_DOM_CUBE); break;
00332   case WAV_CON_CUBE: write_file_class(&outcube, FLA_CON_CUBE); break;
00333   case SUPER_CLASS:  write_file_class(&outcube, SUPER_CLASS);  break;
00334   }
00335 
00336   close_tiger_frame(&incube);
00337   close_tiger_frame(&fcube);
00338   close_tiger_frame(&outcube);
00339 
00340   exit_session(OK);
00341   return(OK);
00342 }

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