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

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