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

focus/source/focus_spectro.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00016 /* =========================================================== */
00017 
00018 #include <IFU_io.h>
00019 #include <IFU_math.h>
00020 #include <snifs.h>
00021 #include <focus.h>
00022 #include "../../extract/incl/extract.h" /* Probably not the right way... */
00023 #include <gsl/gsl_statistics.h>
00024 
00025 #ifdef HAVE_LIBDISLIN
00026 #include <dislin.h>
00027 #endif
00028 
00029 /* === Doxygen Comment ======================================= */
00080 /* =========================================================== */
00081 
00082 int main(int argc, char **argv)
00083 {
00084   int const NARCMAX=20;        /* Maximum number of arcs in input catalogue */
00085   int const ARCWIDTH=128; /* Full-width of arc raster used for correlations */
00086 
00087   IMAGE2D arcframe;
00088   char **argval, **arglabel;
00089   char *catname, *devname, arcname[lg_name+1];
00090   int CAT, MISSSPECFOC=0, SEQ;
00091   int i,j;
00092   int status, narc, foc_start, foc_step;
00093   float specfoc, minfoc, maxfoc, range, ftmp;
00094   double **arctab;
00095   double parpeak[6],eparpeak[6];
00096   double peakx[NARCMAX], *peaky, *peakdy, peakf[NARCMAX];
00097   double peakI[NARCMAX],peakS[NARCMAX],peakSx[NARCMAX],peakSy[NARCMAX];
00098   double peakdI[NARCMAX],peakdS[NARCMAX],peakdSx[NARCMAX],peakdSy[NARCMAX];
00099   double par[3],epar[3];
00100   double msig,emsig,chi2,gof,xnorm;
00101 #ifdef HAVE_LIBDISLIN
00102   Plot plot;
00103   char title[4*lg_ident+1];
00104   float *x,*y,*dy;
00105 #endif
00106 
00107   set_purpose("Compute spectrograph focus from a catalogue of arc exposures");
00108   set_arglist("-in|cat focus.cat -dev ASCII -seq NULL -range 3.");
00109 
00110   init_snifs("$Name:  $","$Revision: 1.22 $");
00111   init_session(argv,argc,&arglabel,&argval);
00112 
00113   if (DEBUG) {
00114     print_msg("$Id: focus_spectro.c,v 1.22 2005/09/21 09:14:27 ycopin Exp $");
00115     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00116   }
00117 
00118   /* ===== Parameters ============================== */
00119 
00120   catname = argval[0];                       /* Input file/catalogue */
00121 
00122   if ((CAT = (strstr(catname,".cat") != NULL))) {
00123     print_msg("o Arc catalog: %s",catname);
00124     RD_catalog(catname, arcname);
00125   } 
00126   else {
00127     strcpy(arcname, catname);
00128     print_msg("o Arc name: %s",arcname);
00129   }
00130 
00131   /* Compare devname=argval[1] to "ASCII" */
00132   if (!strcasecmp((devname = argval[1]),"ASCII") ) {
00133     print_msg("o ASCII output");
00134     *devname = '\0';                         /* Empty string = ASCII */
00135   }
00136   else {
00137 #ifdef HAVE_LIBDISLIN
00138     print_msg("o Output DISLIN device: %s",devname);
00139 #else
00140     print_error("DISLIN output not supported in current executable");
00141     exit_session(ERR_BAD_PARAM);
00142 #endif
00143   }
00144 
00145   if ( (SEQ = is_set(argval[2])) ) {         /* Focus sequence */
00146     if (sscanf(argval[2],"%d,%d",&foc_start,&foc_step) != 2) {
00147       print_error("Invalid option -seq start,step (has '%s')",argval[2]);
00148       exit_session(ERR_BAD_PARAM);
00149     }
00150     else print_msg("o Default linear focus sequence: start = %d, step = %d",
00151                    foc_start, foc_step);
00152   }
00153 
00154   get_argval(3,"%f",&range);                 /* Focus-range quality level */
00155 
00156   /* ===== Loop on catalog (if needed) ============================== */
00157 
00158   alloc2d(&arctab,ARCWIDTH,ARCWIDTH,DOUBLE);
00159   narc = 0;                                  /* Effective nb of arcs */
00160 
00161   do {                                       /* Loop over arc frames */
00162 
00163     /* Open image */
00168     if ((status = open_image(&arcframe,arcname,"I",-1,"arc"))) 
00169       continue;
00170 
00171     /* Extract central raster */
00172     if (extract_arc(&arcframe, ARCWIDTH, arctab)) { 
00173       close_frame(&arcframe);
00174       continue;
00175     }
00176 
00177     /* Read spectrograph position */
00178     specfoc = read_focus(&arcframe);
00179 
00180     if (specfoc < 0) {                       /* Frame number */
00181       MISSSPECFOC = 1;
00182       specfoc = ABS(specfoc);
00183     } 
00184     else {                                   /* Focus position */
00185       if (SEQ && !narc) {                    /* Warn at first frame */
00186         print_warning("Cannot force focus position, option -seq discarded");
00187         SEQ = 0;
00188       }
00189       
00190       if (MISSSPECFOC) {
00191         print_error("Cannot mix focus position and frame number");
00192         close_frame(&arcframe);
00193         exit_session(UNKNOWN);
00194       }
00195     }
00196 
00197     /* Spectrograph position normalization by the first spectrograph position */
00198     if (!narc) {                             /* First spectrograph position */
00199       if (specfoc) {                         /* Non-zero position */
00200         xnorm = specfoc;
00201         if (DEBUG) print_msg("Spectrograph position normalization: %.2f",xnorm);
00202       }
00203       else {                                 /* Zero position */
00204         print_warning("Cannot properly normalize spectrograph position");
00205         xnorm = 1;
00206       }
00207     }
00208     peakx[narc] = specfoc/xnorm;             /* Normalization */
00209 
00210     /* Close arc frame */
00211     close_frame(&arcframe);
00212 
00213     /* Compute and fit auto-correlation peak */
00214     if (DEBUG) print_msg("Computing and fitting arc "
00215                          "auto-correlation peak (%d)",narc+1);
00216 
00217     if (Fit_Xpeak(arctab,arctab,ARCWIDTH,ARCWIDTH,10,ARCWIDTH/4,
00218                   parpeak,eparpeak)) continue;
00219 
00220     if (DEBUG) print_parpeak(parpeak, eparpeak);
00221 
00222     /* Convert correlation quantities to intrinsic quantities */
00223     comp_mean_disp(parpeak, eparpeak, &msig, &emsig);
00224 
00225     /* Store results in growing arrays */
00226     peakI[narc]  = parpeak[2]; peakdI[narc]  = eparpeak[2];
00227     peakSx[narc] = parpeak[3]; peakdSx[narc] = eparpeak[3];
00228     peakSy[narc] = parpeak[4]; peakdSy[narc] = eparpeak[4];
00229     peakS[narc]  = msig;       peakdS[narc]  = emsig;
00230 
00231     narc++;
00232 
00233     if (narc == NARCMAX) {
00234       print_warning("Cannot handle more than %d arcs, keep going",NARCMAX);
00235       break;
00236     }
00237 
00238   } while (CAT && RD_catalog(catname, arcname)); /* Loop over arc frames */
00239 
00240   free2d(&arctab,DOUBLE);
00241 
00242   /* ===== Focus Output ============================== */
00243 
00244   print_msg(" #   Focus    Intens     SigX     SigY    <Sig>   d<Sig>");
00245   for (i=0; i<narc; i++) 
00246     print_msg("%2d %7.0f %9.3g %8.3f %8.3f %8.3f %8.2g",
00247               i+1,peakx[i]*xnorm,
00248               peakI[i],peakSx[i],peakSy[i],peakS[i],peakdS[i]);
00249 
00250   /* Which quantity to be *minimized*? */
00251   peaky  = peakS; peakdy = peakdS;
00252   print_msg("Focus indicator: <Sig> = RMS(SigX,SigY)");
00253 
00254   if (!CAT) {
00255     /* Minimal output */
00256     if (!VERBOSE) {
00257       printf("%f\n",peaky[0]);
00258       fflush(stdout);
00259     }
00260     exit_session(OK);
00261   }  
00262 
00263 #ifdef HAVE_LIBDISLIN
00264   if (strlen(devname) && CAT) {
00265     plot_initialize(&plot,devname,NULL,NULL,0,0);
00266     
00267     /* Conversion from double to float */
00268     x = (float *)malloc(narc*sizeof(float));
00269     y = (float *)malloc(narc*sizeof(float));
00270     dy= (float *)malloc(narc*sizeof(float));
00271     for (i=0; i<narc; i++) {
00272       x[i] = peakx[i]*xnorm; 
00273       if (SEQ) x[i] = foc_start + (x[i]-1)*foc_step;
00274       y[i]  = peaky[i]; 
00275       dy[i] = peakdy[i];
00276     }
00277 
00278     plot_set_limits(x,y,narc,&(plot.Limits));
00279     plot_increase_limits(&(plot.Limits));
00280 
00281     plot_set_device(plot.devname);           /* Output device or filename */
00282     disini();
00283     if (DEBUG) errmod("ALL","ON");
00284     else errmod("ALL","OFF");
00285     complx(); texmod("ON");                  /* TeX support */
00286     labdig(2,"Y");                   /* 2 digits on Y-labels (instead of 1) */
00287     if (VERBOSE) paghdr("$Name:  $", "",1,0); /* Technical header */
00288     plot_axes(&plot,"Focus","<\\sigma >");
00289     sprintf(title,"Spectrograph Focus\nCatalogue: %s\n",catname);
00290     chncrv("BOTH");                 /* Automatic change of curve attributes */
00291     plot_array_err(&(plot.Symbol),x,y,narc,NULL,dy);
00292   }
00293 #endif
00294 
00295   /* ===== Focus fit ============================== */
00296 
00297   /* Initial guesses */
00298   if ( peaky[1] < peaky[0] ) {               /* Concave parabola */
00299     i = gsl_stats_min_index(peaky,1,narc);
00300     j = gsl_stats_max_index(peaky,1,narc);
00301   }
00302   else {                                     /* Convex parabola */
00303     i = gsl_stats_max_index(peaky,1,narc);
00304     j = gsl_stats_min_index(peaky,1,narc);
00305   }
00306   if (peakx[i] == peakx[j]) {
00307     print_error("All the frames seem to have the same focus, aborting");
00308     exit_session(UNKNOWN);
00309   }
00310   par[1] = peakx[i];                         /* Extremum position */
00311   par[2] = peaky[i];                         /* Extremum height */
00312   par[0] = (peaky[j]-par[2])/SQ(peakx[j] - par[1]); /* MinMax parabola */
00313 
00314   if (DEBUG) {
00315     print_msg("Focus guess: %g px at focus = %f [curv=%g]",
00316               par[2],par[1],par[0]);
00317 #ifdef HAVE_LIBDISLIN
00318     if (strlen(devname) && CAT) {
00319       for (i=0; i<narc; i++) y[i] = par[0]*SQ(x[i]/xnorm-par[1]) + par[2];
00320       curve(x,y,narc);                       /* Estimated parabola */
00321     }
00322 #endif
00323   }
00324 
00325   epar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00326 
00327   status = nllsqfit_bnd(peakx, peaky, peakdy, narc, par, NULL, NULL, 3,
00328                         *nllsqfit_quad, peakf, NULL, epar, &chi2,&gof);
00329 
00330   par[0] /= SQ(xnorm); epar[0] /= SQ(xnorm); /* Renormalize */
00331   par[1] *= xnorm;     epar[1] *= xnorm;
00332 
00333   ftmp = sqrt( MAX(0, range/100.*par[2]/par[0]) ); /* Avoid NaN */
00334   minfoc = par[1] - ftmp;       /* Focus range */
00335   maxfoc = par[1] + ftmp;
00336 
00337   if (SEQ) {
00338     print_msg("Converting frame number to focus position (%d,%d)",
00339               foc_start,foc_step);
00340     par[1] = foc_start + (par[1]-1)*foc_step;
00341     epar[1]=                epar[1]*foc_step;
00342     minfoc = foc_start + (minfoc-1)*foc_step;
00343     maxfoc = foc_start + (maxfoc-1)*foc_step;
00344   }
00345 
00346   print_msg("Focus fit: status=%d, Chi2=%g, GoF=%g",status,chi2,gof);
00347   print_msg("Focus position: %.2f +- %.2f",par[1],epar[1]);
00348   print_msg("Best focus: %.2f +- %.2f px",par[2],epar[2]);
00349   print_msg("Focus range [%.0f%%]: %.2f - %.2f",range,minfoc,maxfoc);
00350   if (DEBUG) print_msg("Focus curvature: %g +- %g",par[0],epar[0]);
00351 
00352   /* Minimal output */
00353   if (!VERBOSE) {
00354     printf("%f\n",par[1]);
00355     fflush(stdout);
00356   }
00357 
00358 #ifdef HAVE_LIBDISLIN
00359   if (strlen(devname) && CAT) {
00360     float x0=par[1],y0=par[2],dx0=epar[1],dy0=epar[2];
00361 
00362     for (i=0; i<narc; i++) y[i] = peakf[i]; /* double -> float conversion */
00363 
00364     curve(x,y,narc);                         /* Adjusted parabola */
00365 
00366     plot.Symbol.shape = 21;                  /* Best focus */
00367     plot_array_err(&(plot.Symbol),&x0,&y0,1,&dx0,&dy0);
00368 
00369     color("BLUE");
00370     rline(minfoc,par[2]*(1+range/100.),      /* Focus range */
00371           maxfoc,par[2]*(1+range/100.));
00372     color("FORE");
00373 
00374     free(x); free(y); free(dy);
00375 
00376     sprintf(title,"%sBest focus: %.2f \\pm  %.2f (%.2f \\pm  %.2f)\n",
00377             title,x0,dx0,y0,dy0);
00378     plot_title(title);
00379 
00380     if (strcmp(devname,"CONS")) print_msg("Output file: %s",getfil());
00381 
00382     disfin();
00383   }
00384 #endif
00385 
00386   /* Test if optimal focus is inside the probbed range */
00387   if (par[1]<minfoc || par[1]>maxfoc) {
00388     print_error("Focus position (%.2f) outside %.0f%%-range [%.0f,%.0f]",
00389                 par[1],range,minfoc,maxfoc);
00390     status = 5;                              /* Unused by nllsqfit_bnd */
00391   }
00392 
00393   /* Do not crash if status = 1 
00394      (Optimal solution found, but requested accuracy not achieved) */
00395   if (status == 1) status = 0;
00396 
00397   exit_session(status);
00398   return(status);
00399 }

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