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

extract/source/find_max.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00022 /* =========================================================== */
00023 
00024 #include <IFU_io.h>
00025 #include <IFU_math.h>
00026 #include <snifs.h>
00027 #include <extract.h>
00028 
00029 #include <gsl/gsl_sf_erf.h>
00030 
00032 #define DBG_PREFIX "dbg_fmx"
00033 
00034 /* Global variables describing the geometric PSF */
00035 int glnggeo;                                 /* Nb of gaussian components */
00036 double *glIgeo, *glXgeo, *glSgeo;
00037 
00038 /* Global variables describing the global PSF */
00039 int glngpsf;                                 /* Nb of gaussian components */
00040 
00041 /* Global variables describing the core PSF = geom. PSF (x) global PSF */
00042 double **glIpsf, **glSpsf;           /* Global variables set by set_corePSF */
00043 
00044 /* Global variables describing the local PSF */
00045 int glngloc;                                 /* Nb of gaussian components */
00046 
00047 /* === Doxygen Comment ======================================= */
00069 /* =========================================================== */
00070 
00071 int main(int argc, char **argv)
00072 {
00073   Channel channel, channel2;
00074 
00075   IMAGE2D contima;
00076   TABLE mask;
00077   SPECTRUM fit, med;
00078   Maxima_Set maxset;
00079 
00080   char **argval, **arglabel, text[132];
00081   char tmpname[lg_name+1];
00082   int i,j,k,l, jstart, jend, singleline, status;
00083   int **imax, **imin, *nmax, *iwork, npts;
00084   int imedline, nmedline, i0, i1, nlmin, nlmax, Ntmin, Ntmax;
00085   int fullblock, halfblock, rejpeak;
00086   int fitpeak, npeak, totpeak;
00087   int curpeak, npack, packpb, packpe;
00088   int nggeo, ngpsf, ngloc;
00089 
00090   float  y, xstart;
00091   double scale, *xprof, *xfit;
00092   double *intens, **median_line, threshold, rms;
00093   double *psf_param, *geo_param, *packpar, *dpackpar;
00094   double *fmin, *fmax;
00095   double mean,sigma,meanmed,sigmed;
00096 
00097   /* ##### INTRODUCTION ############################## */
00098 
00099   set_arglist("-in none -max max -linestep 5 -mask mask -singleline 0");
00100   set_purpose("Search for maxima of X-dispersion profile in continuum frame");
00101 
00102   init_snifs("$Name:  $");
00103   init_session(argv,argc,&arglabel,&argval);
00104 
00105   if (DEBUG) {
00106     print_msg("$Id: find_max.c,v 1.16 2004/10/06 16:56:32 ycopin Exp $");
00107     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00108   }
00109 
00110   /* ===== Parameters ============================== */
00111 
00112   get_argval(2,"%d",&fullblock); /* Nb of lines to average */
00113 
00114   print_msg("o Input continuum frame: %s",argval[0]);
00115   print_msg("o Output max file: %s",argval[1]);
00116   print_msg("o Nb of lines to average in analysis: %d",fullblock);
00117   print_msg("o Input mask: %s",argval[3]);
00118 
00119   get_argval(4,"%d",&singleline); /* Single line fit */
00120   if (singleline) print_warning("Processing line %d only", singleline);
00121 
00122   /* ===== Mask ============================== */
00123 
00124   if ((status = open_mask(argval[3],&mask,"I"))) exit_session(status);
00125   if (!(channel = read_channel(&mask))) exit_session(channel);
00126 
00127   /* ===== PSF models ============================== */
00128 
00129   print_msg("Reading PSF models from mask descriptors");
00130 
00131   /* ----- Geometric PSF ------------------------------ */
00132 
00133   if (RD_desc(&mask,MPUP_NGGEO,INT,1,&nggeo) < 0) {
00134     print_error("Unable to read descriptor %s",MPUP_NGGEO);
00135     exit_session(ERR_NODESC);
00136   }
00137   geo_param = (double *)malloc(3*nggeo*sizeof(double));
00138   if (RD_desc(&mask,MPUP_GEOPAR,DOUBLE,3*nggeo,geo_param) < 0) {
00139     print_error("Unable to read descriptor %s",MPUP_GEOPAR);
00140     exit_session(ERR_NODESC);
00141   }
00142   print_msg("   Geometric PSF model: %d component(s)",nggeo);
00143 
00144   glnggeo = nggeo;
00145   glIgeo = (double *)malloc(nggeo*sizeof(double));
00146   glXgeo = (double *)malloc(nggeo*sizeof(double));
00147   glSgeo = (double *)malloc(nggeo*sizeof(double));  
00148   /* GEOPAR already include the scale parameter, and is properly normed */
00149   for (i=0; i<nggeo; i++) {
00150     glIgeo[i] = geo_param[3*i  ];
00151     glXgeo[i] = geo_param[3*i+1];
00152     glSgeo[i] = geo_param[3*i+2];
00153     if (DEBUG) print_msg("\tGaussian %d: I = %f, X = %f, S = %f", 
00154                          i+1,geo_param[3*i],geo_param[3*i+1],geo_param[3*i+2]);
00155   }
00156   free(geo_param);
00157 
00158   /* ----- Global PSF ------------------------------ */
00159 
00160   if (RD_desc(&mask,MPUP_NGPSF,INT,1,&ngpsf) < 0) {
00161     print_error("Unable to read descriptor %s",MPUP_NGPSF);
00162     exit_session(ERR_NODESC);
00163   }
00164   psf_param = (double *)malloc((1+2*ngpsf)*sizeof(double));
00165   if (RD_desc(&mask,MPUP_PSFPAR,DOUBLE,(1+2*ngpsf),psf_param) < 0) {
00166     print_error("Unable to read descriptor %s",MPUP_PSFPAR);
00167     exit_session(ERR_NODESC);
00168   }
00169   scale = psf_param[0];
00170   print_msg("   Geometric pupil scale: %8.4f pix", scale);
00171   print_msg("   Global PSF model: %d component(s)",ngpsf);
00172 
00173   if (DEBUG) 
00174     for (i=0; i<ngpsf; i++) 
00175       print_msg("\tUnconvolved gaussian %d: I = %f, S = %f", 
00176                 i+1,psf_param[2*i+1],psf_param[2*i+2]);
00177 
00178   /* ----- Local PSF ------------------------------ */
00179 
00180   if (RD_desc(&mask,MPUP_NGLOC,INT,1,&ngloc) < 0) {
00181     print_error("Unable to read descriptor %s",MPUP_NGLOC);
00182     exit_session(ERR_NODESC);
00183   }
00184   print_msg("   Local PSF model: %d component(s)",ABS(ngloc));
00185 
00186   close_table(&mask);
00187 
00188   /* We do not handle local PSF with more than 1 gaussian... */
00189   if (ABS(ngloc) > 1) {
00190     print_error("Local PSF with more than *one* gaussian not implemented yet!");
00191     exit_session(ERR_NOIMPL);
00192   }
00193 
00194   /* ----- Core PSF ------------------------------ */
00195 
00196   if (ngloc > 0) {
00197     print_msg("   Core PSF = Geometric PSF (x) Global PSF");
00198   }
00199   else {
00200     ngloc = -ngloc;
00201     ngpsf = 0;
00202     print_msg("   Core PSF = Geometric PSF (global PSF *discarded*)");
00203   }
00204 
00205   if (ngloc > 1) {
00206     print_error("Multi-gaussian local PSF model not implemented yet");
00207     exit_session(ERR_NOIMPL);
00208   }
00209 
00210   glngpsf = ngpsf;
00211   alloc2d(&glIpsf,nggeo,MAX(ngpsf,1),DOUBLE);
00212   alloc2d(&glSpsf,nggeo,MAX(ngpsf,1),DOUBLE);
00213 
00214   /* If ngpsf=0, psfpar is not needed, since core PSF = geom. PSF then */
00215   set_corePSF(psf_param);
00216 
00217   free(psf_param);
00218 
00219   /* ===== Continuum frame ============================== */
00220 
00221   if ((status = open_image(&contima,argval[0],"I",-1,"continuum")))
00222     exit_session(status);
00223 
00224   /* Check channel */
00225   if ((channel2 = read_channel(&contima)) != channel) {
00226     print_error("mask channel (%s) and input frame channel (%s) are different",
00227                 channel_names[channel],channel_names[channel2]);
00228     exit_session(ERR_BAD_PARAM);
00229   }
00230 
00231   /* ----- Frame statistics ------------------------------ */
00232 
00233   //  frame_stats(&contima,&mean,&sigma,&meanmed,&sigmed);
00234   //  print_msg("   Mean=%.5f (%.5f), Median=%.5f (%.5f)",
00235   //            mean,sigma,meanmed,sigmed);
00236 
00237   /* ##### COMPUTATIONS ############################## */
00238 
00239   /* ===== Search for maxima ============================== */
00240 
00241   intens = (double *)malloc(fullblock*sizeof(double));
00242   iwork  =    (int *)malloc(fullblock*sizeof(int));
00243 
00244   halfblock = (int)(fullblock/2.);             /* half block size */
00245   jstart = halfblock;                          /* minimal central line (incl)*/
00246   jend   = contima.ny - fullblock + halfblock; /* maximal central line (excl)*/
00247   jend   = MAX(jstart+1,jend);
00248 
00249   //#ifdef SINGLELINE_THRESHOLD
00250   if (singleline) {
00251     /* Limit the threshold computation to the current block. Much
00252        faster, but may change the value used (not real DEBUG) */
00253     jstart = MAX(singleline,  jstart);
00254     jend   = MIN(singleline+1,jend);
00255     print_warning("Computing threshold from single line block @[%d-%d[",
00256                   jstart-halfblock,jstart-halfblock+fullblock);
00257   }
00258   //#endif
00259 
00260   /* Total nb of central lines=blocks to process */
00261   nmedline = (int)((jend-1 - jstart)/fullblock) + 1; 
00262 
00263   /* Median line */
00264   alloc2d(&median_line,nmedline,contima.nx,DOUBLE);
00265   /* Index of min and max in median line */
00266   alloc2d(&imin,nmedline,contima.nx/2.,INT);
00267   alloc2d(&imax,nmedline,contima.nx/2.,INT);
00268   /* Flux of min and max in median line */
00269   fmin = (double *)malloc(nmedline*contima.nx/2.*sizeof(double));
00270   fmax = (double *)malloc(nmedline*contima.nx/2.*sizeof(double));
00271   /* Number of max in median line */
00272   nmax = (int *)malloc(nmedline*sizeof(int));
00273 
00274   Ntmax = 0; /* Total number of maxima & minima */
00275   Ntmin = 0; 
00276 
00277   reset_print_progress();
00278   for (j=jstart; j<jend; j+=fullblock) { /* loop over continuum rows */
00279 
00280     /* DEF: imedline: median line index (0..nmedline-1) */
00281     imedline = (int)((j - jstart)/fullblock);
00282 
00283     if (!DEBUG) print_progress("Thresholding",(imedline+1.)/nmedline*100.,5.);
00284 
00285     /* ----- Median of fullblock rows ------------------------------ */
00286 
00287     for (i=0; i<contima.nx; i++) { /* Loop along the central line */
00288       for (l=0; l<fullblock; l++)  /* Loop on the lines of the block */
00289         intens[l] = RD_frame(&contima,i,j-halfblock+l);
00290       median_line[imedline][i] = median(intens,fullblock,iwork);
00291     }
00292 
00293     /* ----- Extrema search ------------------------------ */
00294 
00295     nlmax = search_extrema(median_line[imedline], contima.nx, 
00296                            imin[imedline], fmin+Ntmin, 
00297                            imax[imedline], fmax+Ntmax, &nlmin);
00298     
00299     nmax[imedline] = nlmax; /* number of max detected in current line */
00300 
00301     /* We might miss the last minimum, not needed anyway */
00302 
00303     if (DEBUG) {
00304       mean_sig(fmax+Ntmax,nlmax,&mean,&sigma);
00305       mean_sig(fmin+Ntmin,nlmin,&meanmed,&sigmed); /* temp. var. */
00306       print_msg("Block #%d @[%d-%d[: %d max (%.3g +- %.3g), "
00307                 "%d min (%.3g +- %.3g)",
00308                 imedline+1,j-halfblock,j-halfblock+fullblock,
00309                 nlmax,mean,sigma,nlmin,meanmed,sigmed);
00310     }
00311 
00312     /* Update */
00313     Ntmax += nlmax;
00314     Ntmin += nlmin;
00315 
00316   } /* continuum rows */
00317 
00318   print_msg("Total nb of maxima detected: %d (%.1f max/line)", 
00319             Ntmax,1.*Ntmax/nmedline);
00320 
00321   free(intens); free(iwork);
00322   close_frame(&contima);
00323 
00324   /* ----- Threshold ------------------------------ */
00325 
00326   /* compute statistics on minima to discard maxima at noise level */
00327   if (Ntmin == 0) {
00328     print_error("ARGH! No minimum found in frame");
00329     exit_session(ERR_BAD_TYPE);
00330   }
00331   mean_sig(fmin,Ntmin,&mean,&sigma); /* overall stats on minima */
00332   if (DEBUG) print_msg("Total nb of minima detected: %d (%g +- %g)", 
00333                        Ntmin,mean,sigma);
00334   threshold = mean + FMX_SIGMA_CLIP*sigma;
00335   print_msg("Computed Threshold (%.1f-sigma): %g", FMX_SIGMA_CLIP,threshold);
00336 
00337   print_msg("Neighborhood limit: %.1f x %.1f = %.1f pix",
00338             FMX_TOOFAR,SNIFS_INTERSP,FMX_TOOFAR*SNIFS_INTERSP);
00339 
00340   free(fmin); free(fmax);
00341 
00342   /* ===== Cross-dispersion profiles fitting ============================== */
00343 
00344   /* ----- Max structure ------------------------------ */
00345 
00346   if (!singleline) {
00347     if (alloc_TIGER_max(&maxset,nmedline,(int)(contima.nx/2)) != 0) {
00348       print_error("Cannot allocate max structure");
00349       exit_session(ERR_ALLOC);
00350     }
00351     maxset.nb_ycoords = nmedline;
00352   }
00353   else {
00354     imedline = (int)(singleline - jstart)/fullblock;
00355 
00356     print_warning("Working on single line %d (block #%d @[%d-%d[): %d max",
00357                   singleline,imedline+1,
00358                   singleline-halfblock,singleline-halfblock+fullblock,
00359                   nmax[imedline]);
00360 
00361     /* Storage of median line and fit */
00362     sprintf(tmpname,"%s_med",DBG_PREFIX);
00363     print_msg("DEBUG: Creating %s with cross-dispersion profile",tmpname);
00364     disable_erase_flag();
00365     create_spec(&med,tmpname,contima.nx,contima.startx,
00366                 contima.stepx,DOUBLE,"curve to fit","pixels");
00367     for (i=0; i<med.npts; i++) WR_spec(&med,i,median_line[imedline][i]);
00368     write_file_class(&med, DEBUG_FILE);
00369     close_spec(&med);
00370 
00371     sprintf(tmpname,"%s_fit",DBG_PREFIX);
00372     print_msg("DEBUG: Creating %s with cross-dispersion fit",tmpname);
00373     create_spec(&fit,tmpname,contima.nx,contima.startx,
00374                 contima.stepx,DOUBLE,"fitted curve","pixels");
00375 
00376     /* Get ready for loop on lines */
00377     jstart = singleline;
00378     jend   = singleline+1;
00379   }
00380 
00381   /* ----- Loop on lines ------------------------------ */
00382 
00383   totpeak = 0; /* nb of peaks to be adjusted */
00384   rejpeak = 0; /* nb of rejected peaks */
00385 
00386   reset_print_progress();
00387   for (j=jstart; j<jend; j+=fullblock) { /* Loop on central lines */
00388 
00389     /* DEF: imedline: median line index */
00390     imedline = (int)((j - jstart)/fullblock);
00391 
00392     coord_frame(&contima,0,j,&xstart,&y);
00393 
00394     if (!singleline && !DEBUG) {
00395       sprintf(text,"Profile Fitting: line %d (y=%.1f), block #%d @[%d-%d[",
00396               j,y,imedline+1,j-halfblock,j-halfblock+fullblock);
00397       print_progress(text, 100*(imedline + 1)/nmedline, 1.0);
00398     }
00399     else {
00400       print_msg("Working on line %d (y=%.1f, block #%d @[%d-%d[): %d max",
00401                 j,y,imedline+1,j-halfblock,j-halfblock+fullblock,
00402                 nmax[imedline]);
00403     }
00404 
00405     curpeak = 0; /* index of current peak in current line */
00406     fitpeak = 0; /* nb of adjusted peak in current line */
00407     npack   = 0; /* nb of packs in current line */
00408 
00409     while (curpeak < nmax[imedline]) { /* Loop on the peaks of the line */
00410 
00411       if (singleline && !DEBUG) {
00412         sprintf(text,"Profile Fitting: peak %d/%d",curpeak,nmax[imedline]);
00413         print_progress(text,(curpeak+1.)/nmax[imedline]*100.,1.);
00414       }
00415       
00416       npeak = search_pack(median_line[imedline], &curpeak, 
00417                           imax[imedline], nmax[imedline], 
00418                           threshold, &packpb, &packpe);
00419 
00420       if (!npeak) break; /* No pack found, stop there */      
00421 
00422       totpeak += npeak;
00423       npack ++; 
00424 
00425       if (singleline && DEBUG) {
00426         print_msg("   Fitting pack #%d with %d peaks (%d to go): "
00427                   "peaks %d (@%d) to %d (@%d)",
00428                   npack,npeak,nmax[imedline]-packpe-1,
00429                   packpb+1,imax[imedline][packpb],
00430                   packpe+1,imax[imedline][packpe]);
00431       }
00432 
00433       /* Limits of the fit: we don't have the minimum after the last max */
00434       i0 = imin[imedline][packpb];
00435       if (packpe<nmax[imedline]-1) i1 = imin[imedline][packpe+1];
00436       else                         i1 = contima.nx - 1;
00437       npts = i1-i0+1;
00438 
00439       /* Storage */
00440       xprof   = (double *)malloc(   npts*sizeof(double));
00441       xfit    = (double *)malloc(   npts*sizeof(double));
00442       packpar = (double *)malloc(3*npeak*sizeof(double));
00443       dpackpar= (double *)malloc(3*npeak*sizeof(double));
00444       for (i=0; i<npts; i++) xprof[i] = median_line[imedline][i0+i];
00445 
00446       /* Initial guesses */
00447       for (k=0; k<npeak; k++) {
00448         /* For position: center @ max */
00449         packpar[3*k  ] = imax[imedline][packpb+k] - i0;
00450         /* For width: HARD-CODED!!! */
00451         packpar[3*k+2] = 0.1;
00452         /* For intensity */
00453         packpar[3*k+1] = xprof[((int)packpar[3*k])]/(packpar[3*k+2]*M_SQRT2PI);
00454 
00455 #ifdef MEGAVERBOSE
00456         if (singleline) {
00457           print_msg("         Peak %2dE: X=%.1f\tI=%7.1e\tS=%.2f",
00458                     k+1,packpar[3*k]+i0 + contima.startx,
00459                     packpar[3*k+1],packpar[3*k+2]);
00460         }
00461 #endif
00462       }
00463 
00464       rms = Fit_pack(xprof,npts,npeak,packpar,dpackpar,xfit);
00465 
00466 #ifdef MEGAVERBOSE
00467       if (singleline) {
00468         if (rms>0) {
00469           print_msg("         RMS = %.3g",rms);
00470           for (k=0; k<npeak; k++) {
00471             print_msg("         Peak %2d:  X=%.1f\tI=%7.1e\tS=%.2f",
00472                       k+1,packpar[3*k]+i0 + contima.startx,
00473                       packpar[3*k+1],packpar[3*k+2]);
00474           }
00475         }
00476         else print_msg("         Status %d", status);
00477       }
00478 #endif
00479 
00480       if (rms>0) { /* "Successful" fitting of the npeak peak in the pack */
00481         if (singleline) { /* Storage in fit spectrum */
00482           for (i=0; i<npts; i++) WR_spec(&fit,i0+i,xfit[i]);
00483         }
00484         else {            /* Storage in maximum structure */
00485           for (k=0; k<npeak; k++) {
00486             maxset.line[imedline].maxima[fitpeak+k].xcoord  = 
00487               packpar[3*k]+i0 + contima.startx; /* we switch back to coord. */
00488             maxset.line[imedline].maxima[fitpeak+k].intens  = packpar[3*k+1];
00489             maxset.line[imedline].maxima[fitpeak+k].sigma[0]= packpar[3*k+2];
00490             /* Let's store the error on x in alpha */
00491             maxset.line[imedline].maxima[fitpeak+k].alpha   = dpackpar[3*k];
00492           }
00493         }
00494         fitpeak += npeak; /* Update nb of successfully adjusted peaks */
00495       } else      /* Unsuccessful fit */
00496         rejpeak += npeak; /* Update nb of un-successfully adjusted peaks */
00497 
00498       free(xprof); free(xfit);
00499       free(packpar); free(dpackpar);
00500 
00501     } /* Loop on the peaks of the line */
00502     
00503     if (singleline) {
00504       print_msg("Line %d: %d peaks detected, %d/%d fitted in %d packs", 
00505                 j,nmax[imedline],fitpeak,totpeak,npack);
00506     }
00507     else {
00508       maxset.line[imedline].ycoord = y;
00509       maxset.line[imedline].nb_max = fitpeak;
00510     }
00511 
00512   } /* Loop on central lines */
00513 
00514   print_msg("Fraction of rejected peaks: %6.2f%% (%d/%d)",
00515             100.*rejpeak/totpeak,rejpeak,totpeak);
00516 
00517   if (singleline) {
00518     write_file_class(&fit, DEBUG_FILE);
00519     close_spec(&fit);
00520   }
00521   else {
00522     print_msg("Saving max structure to %s", argval[1]);
00523     status = save_TIGER_max(&maxset,argval[1]);
00524   }
00525 
00526   free2d(&glIpsf,DOUBLE); free2d(&glSpsf,DOUBLE);
00527   free(glIgeo); free(glXgeo); free(glSgeo);
00528 
00529   free2d(&imin,INT); free2d(&imax,INT);
00530 
00531   exit_session(OK);
00532   return(OK);
00533 }
00534 

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