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

calib/source/truncate_cube.c

Go to the documentation of this file.
00001 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00002 ! COPYRIGHT    (c)  1995 CRAL
00003 ! IDENT        truncate_tiger.c
00004 ! LANGUAGE     C
00005 ! AUTHOR       RB
00006 !            
00007 ! PURPOSE      Truncate a tiger cube to a common wavelength
00008 !              All spectra wich do not fit in the interval are rejected
00009 ! COMMENT
00010 !       **** INFO mode (no output file created)
00011 !            truncate_tiger -in file -wave l1,l2 -info [-infotbl tbl]
00012 ! 
00013 !        if -infotbl name option is used, a table is created with
00014 !        a copy of the :No and :A,:D (or :X0,:Y0) columns of the
00015 !        table associated to the input tiger file, and a column :OUT
00016 !        which is set to 1 if the spectrum is rejected. This table
00017 !        can be used to display the spectrum rejected by this command
00018 ! 
00019 !       **** EXECUTE mode (output file created)
00020 !            truncate_tiger -in file -wave l1,l2 -out file
00021 ! 
00022 !     in both case -wave l1,l2 can be replaced by -ref spec and the wavelength
00023 !     range l1,l2 is read from the reference spectrum 
00024 ! 
00025 ! VERSION 0.0  01-09-95 RB: creation
00026 ! VERSION 0.1  15-12-95 RB: nouvelle version IOLIB
00027 ! VERSION 0.2  10-01-96 RB: mise a jour iolib
00028 ! VERSION 0.3  20-12-99 YC: application a un spectre
00029 ! VERSION 1.0  03-06-00 RB: include noise if available
00030 ! VERSION 1.1  25-08-00 RB: add CP_non_std in spec
00031 ! VERSION 1.2  ???
00032 ! VERSION 1.3  14-12-04 YC: close_tiger
00033 ! VERSION 1.4  16-01-05 YC: info, clipping
00034 ---------------------------------------------------------------------*/
00035 
00036 #include <IFU_io.h>
00037 
00038 int get_coordcol_id(TABLE *table, int *colx, int *coly, char *type)
00039 {
00040   set_control_level(NONE);
00041 
00042   *colx = get_col_ref(table, "A");
00043   *coly = get_col_ref(table, "D");
00044   *type = 'A';
00045 
00046   if ((*colx < 0) || (*coly < 0)) {
00047     *colx = get_col_ref(table, "X0");
00048     *coly = get_col_ref(table, "Y0");
00049     *type = 'C';
00050   }
00051 
00052   set_control_level(FATAL);
00053   if ((*colx < 0) || (*coly < 0)) return(ERR_BAD_COL);
00054   return(OK);
00055 }
00056 
00057 int main(int argc, char **argv)
00058 {
00059   TABLE intab, outtab;
00060   TIGERfile incube, outcube;
00061   SPECTRUM inspec, outspec, noise, outnoise;
00062 
00063   char  **argval, **argname;
00064   char *inname, *outname, *tabname;
00065   char separator[1], type;
00066   int col_x, col_y, col_no, n, l, out, nkeep, nrej, NOISE;
00067   int i, k, WAVE, REF, INFO, SPEC, SAVE, NOOUT, no, npts;
00068   int outcol_no, outcol_x, outcol_y, outcol_out;
00069   float x, y;
00070   double lim[2], lbda_min, lbda_max, step;
00071 
00072   /*==================================================================
00073     Input parameters reading 
00074     ==================================================================*/
00075 
00076   set_version("1.4");
00077   set_purpose("Truncate spectra of a datacube to a common wavelength range");
00078   set_arglist("-in none -out NULL -wave NULL -ref NULL -info -infotbl NULL -spec");
00079   init_session(argv, argc, &argname, &argval);
00080 
00081   inname = argval[0];
00082   if (!(NOOUT = is_unset(argval[1]))) outname = argval[1];
00083 
00084   /* Wavelength interval */
00085   WAVE = is_set(argval[2]);
00086   REF  = is_set(argval[3]);
00087   if ((WAVE && REF) || (!WAVE && !REF)) {
00088     print_error("A single option -wave or -ref has to be set");
00089     exit_session(ERR_BAD_PARAM);
00090   }
00091 
00092   if (WAVE && (decode_argval_double(argval[2], lim, separator) != 2 ||
00093                lim[1] < lim[0])) {
00094     print_error("Error in input parameter -wave lbda1,lbda2 (got %s)", 
00095                 argval[2]);
00096     exit_session(ERR_BAD_PARAM);
00097   }
00098 
00099   if (REF) {                                 /* Reference spectrum */
00100     print_msg("Reading wavelength range from spectrum %s", argval[3]);
00101     if (open_spec(&inspec, argval[3], "I")<0) {
00102       print_error("Cannot open reference spectrum %s", argval[3]);
00103       exit_session(ERR_OPEN);
00104     }
00105     lim[0] = inspec.start;
00106     lim[1] = inspec.end;
00107     close_spec(&inspec);
00108   }
00109 
00110   INFO = is_true(argval[4]);                 /* Dry mode */
00111   if ((SAVE = is_set(argval[5])))            /* Info output table */
00112     tabname = argval[5];
00113   SPEC = is_true(argval[6]);                 /* Spectrum flag */
00114 
00115   if (NOOUT && !INFO) {
00116     print_error("Parameter -out filename is required in execution mode");
00117     exit_session(ERR_BAD_PARAM);
00118   }
00119 
00120   /*==================================================================
00121     Initialization
00122     ==================================================================*/
00123 
00124   if (SPEC) {
00125     print_msg("Opening input spectrum %s", inname);
00126     if (open_spec(&inspec, inname, "I")<0) {
00127       print_error("Cannot open input spectrum %s", inname);
00128       exit_session(ERR_OPEN);
00129     }
00130     step = inspec.step;
00131   } 
00132   else {
00133     print_msg("Opening input datacube %s", inname);
00134     if (open_tiger_frame(&incube, inname, "I")<0) {
00135       print_error("Cannot open input datacube '%s'", inname);
00136       exit_session(ERR_OPEN);
00137     }
00138     step = incube.step;
00139     NOISE = exist_noise(&incube);
00140 
00141     if (lim[0] == 0 && lim[1] == 0) {        /* Clipping */
00142       print_msg("Computing clipping range");
00143       if (has_common_bounds(&incube)) {
00144         get_common_param(&incube,&npts,&(lim[0]),&(lim[1]));
00145         print_msg("   Common bounds: %.2f - %.2f (%d pixels)",
00146                   lim[0],lim[1],npts);
00147       }
00148       else {
00149         lim[0] = -MAXFLOAT;
00150         lim[1] =  MAXFLOAT;
00151         for (n=0; n<incube.nbspec; n++) {
00152           get_tiger_spec(&incube, &inspec, NULL, incube.signal[n].specId);
00153           lim[0] = MAX(inspec.start,lim[0]);
00154           lim[1] = MIN(inspec.end,lim[1]);
00155         }
00156         print_msg("   Clipping range: %.2f - %.2f",lim[0],lim[1]);
00157       }
00158     }
00159   }
00160 
00161   /* Calcul des limites en fonction du pas */
00162   lbda_min = (int)(lim[0]/step + 0.5)*step;
00163   lbda_max = (int)(lim[1]/step + 0.5)*step;
00164   npts = ((lbda_max - lbda_min)/step + 1.5);
00165   print_msg("Truncated wavelength range: %f - %f (%d pixels)", 
00166             lbda_min, lbda_max, npts);
00167 
00168   /*==================================================================
00169     Single spectrum case
00170     ==================================================================*/
00171 
00172   if (SPEC) {
00173     if (INFO) {
00174       print_msg("Input spectrum range: %f - %f",inspec.start,inspec.end);
00175       close_spec(&inspec);
00176       exit_session(OK);
00177     }
00178 
00179     if (((float)inspec.start > (float)lbda_min) || 
00180         ((float)inspec.end   < (float)lbda_max)) {
00181       print_error("Cannot truncate spectrum [%f-%f] at "
00182                   "requested wavelength [%f-%f]",
00183                   inspec.start,inspec.end,lbda_min,lbda_max);
00184       close_spec(&inspec);
00185       exit_session(ERR_BAD_SIZE);
00186     }
00187 
00188     print_msg("Creating truncated spectrum %s", outname);
00189     create_spec(&outspec, outname, npts, lbda_min,
00190                 inspec.step, inspec.data_type, inspec.ident, inspec.cunit);
00191     for (i=pixel_spec(&inspec, lbda_min), k=0; k<npts; k++, i++)
00192       WR_spec(&outspec, k, RD_spec(&inspec, i));
00193 
00194     CP_non_std_desc(&inspec, &outspec);
00195     close_spec(&inspec);
00196     close_spec(&outspec);
00197     exit_session(OK);
00198   }
00199 
00200   /*==================================================================
00201     Datacube case
00202     ==================================================================*/
00203 
00204   if (INFO) {
00205     print_msg("INFO mode (no truncation applied)");
00206 
00207     if (SAVE) {                              /* Store info in table */
00208       print_msg("Results saved in table %s, cols No,X,Y,OUT", tabname);
00209 
00210       open_table(&intab, incube.table_name, "I");
00211       col_no = get_col_ref(&intab, LAB_COL_NO);
00212       if (get_coordcol_id(&intab, &col_x, &col_y, &type)) {
00213         print_error("Cannot read spatial coordinates from table %s",intab.name);
00214         exit_session(ERR_BAD_COL);
00215       }
00216 
00217       create_table(&outtab, tabname, intab.allrow, 4, 'I', " ");
00218       outcol_no= create_col(&outtab,LAB_COL_NO,LONG,  'N',"I4","Lens number");
00219       outcol_x = create_col(&outtab, "X",      FLOAT, 'N',"F10.4", "X-coord");
00220       outcol_y = create_col(&outtab, "Y",      FLOAT, 'N',"F10.4", "Y-coord");
00221       outcol_out=create_col(&outtab, "OUT",    LONG,  'N',"I2",  "Discarded");
00222     }
00223 
00224     for (n=0, nrej=0; n<incube.nbspec; n++) {
00225       print_progress("Scanning data cube", 100*(n+1)/incube.nbspec, 10.);
00226       no = incube.signal[n].specId;
00227       if (SAVE) {
00228         l = search_in_col(&intab, col_no, &no);
00229         if (RD_tbl(&intab, l, col_x, &x)) continue;
00230         if (RD_tbl(&intab, l, col_y, &y)) continue;
00231       }
00232 
00233       get_tiger_spec(&incube, &inspec, NULL, no);
00234       out = 0;
00235       if (((float)inspec.start > (float)lbda_min) || 
00236           ((float)inspec.end   < (float)lbda_max)) {
00237         nrej++;
00238         out = 1;
00239       }
00240       if (SAVE) {
00241         WR_tbl(&outtab, no, outcol_no, &no);
00242         WR_tbl(&outtab, no, outcol_x, &x);
00243         WR_tbl(&outtab, no, outcol_y, &y);
00244         WR_tbl(&outtab, no, outcol_out, &out);
00245       }
00246       free_spec_mem(&inspec);
00247     }
00248 
00249     nkeep = n - nrej;
00250     print_msg("Number of spectra: %d (Complete) %d (Discarded) %d (Total)",
00251               nkeep, nrej, n);
00252     if (SAVE) close_table(&outtab);
00253 
00254     exit_session(OK);
00255 
00256   } /* INFO */
00257 
00258   print_msg("Creating truncated datacube %s", outname);
00259   if (create_tiger_frame(&outcube, outname, npts, lbda_min, incube.step,
00260                          incube.data_type,incube.table_name,incube.ident,incube.cunit)) {
00261     print_error("Cannot create datacube %s", outname);
00262     exit_session(ERR_CREAT);
00263   }
00264 
00265   reset_print_progress();
00266   for (n=0, nkeep=0, nrej=0; n<incube.nbspec; n++) {
00267     print_progress("Truncating datacube", 100*(n+1)/incube.nbspec, 1.);
00268     no = incube.signal[n].specId;
00269     if (NOISE) get_tiger_spec(&incube, &inspec, &noise, no);
00270     else       get_tiger_spec(&incube, &inspec, NULL, no);
00271     /* Make the comparison in float to avoid rouding errors */
00272     if (((float)inspec.start > (float)lbda_min) || 
00273         ((float)inspec.end   < (float)lbda_max)) {
00274       if (DEBUG) print_warning("lens #%d discarded (%f-%f)",
00275                                no,inspec.start,inspec.end);
00276       nrej++;
00277       free_spec_mem(&inspec);
00278       continue;
00279     }
00280 
00281     init_new_tiger_spec(&outcube, &outspec, npts, lbda_min);
00282     if (NOISE) init_new_tiger_spec(&outcube, &outnoise, npts, lbda_min);
00283 
00284     for (i=pixel_spec(&inspec, lbda_min), k=0; k<npts; k++, i++) {
00285       WR_spec(&outspec, k, RD_spec(&inspec, i));
00286       if (NOISE) WR_spec(&outnoise, k, RD_spec(&noise, i));
00287     }
00288     if (NOISE) {
00289       put_tiger_spec(&outcube, &outspec, &outnoise, no);
00290       free_spec_mem(&inspec);
00291       free_spec_mem(&noise);
00292     } 
00293     else {
00294       put_tiger_spec(&outcube, &outspec, NULL, no);
00295       free_spec_mem(&inspec);
00296     }
00297     nkeep++;
00298   }
00299   n = nkeep + nrej;
00300   print_msg("Number of spectra: %d (Complete) %d (Discarded) %d (Total)",
00301             nkeep, nrej, n);
00302 
00303   if (nkeep < n/10) 
00304     /* Issue a warning if less than 10% of the spectra went through */
00305     print_warning("Only %.1f%% of the spectra were truncated",nkeep*100./n);
00306   
00307   CP_non_std_desc(&incube, &outcube);
00308   close_tiger_frame(&incube);
00309   close_tiger_frame(&outcube);
00310 
00311   exit_session(OK);
00312   return(OK);
00313 }

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