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

extract/lib/proc_max.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00013 /* =========================================================== */
00014 
00015 #include <IFU_io.h>
00016 #include <IFU_math.h>
00017 #include <snifs.h>
00018 #include <extract.h>
00019 
00020 /* === Doxygen Comment ======================================= */
00035 /* =========================================================== */
00036 
00037 int search_extrema(double line[], int npts, 
00038                    int imin[], double fmin[], int imax[], double fmax[], 
00039                    int *nlmin)
00040 {
00041   int nmax=0, nmin=0, i=0; /* In SAURON, i=1, why? */
00042 
00043   while (i < npts) {
00044     /* Search for a minimum: going down (we might miss the 1st max) */
00045     while (i < npts-1 && line[i] >= line[i+1]) i++;
00046     if (i >= npts-1) break;
00047 
00048     imin[nmin] = i;
00049     fmin[nmin] = line[i];
00050     nmin ++;
00051 
00052     /* Search for the maximum following previous minimum: going up */
00053     while (i < npts-1 && line[i] <= line[i+1]) i++;
00054     if (i >= npts-1) break;
00055     imax[nmax] = i;
00056     fmax[nmax] = line[i];
00057     nmax ++;
00058 
00059     /*
00060       if (DEBUG) 
00061         print_msg("   Ext.%d: Min[%d] = %.3g, Max[%d] = %.3g",
00062                   nmin,imin[nmin-1],fmin[nmin-1],imax[nmax-1],fmax[nmax-1]);
00063     */
00064   }
00065   *nlmin = nmin;
00066   return(nmax);
00067 }
00068 
00069 /* === Doxygen Comment ======================================= */
00082 /* =========================================================== */
00083 
00084 int search_pack(double line[], int *curpeak, int imax[], int nmax, 
00085                 double threshold, int *packpb, int *packpe)
00086 {
00087   int npeak;
00088   
00089   /* Find the 1st valid peak (above threshold) of the pack */
00090   while (*curpeak < nmax && line[imax[*curpeak]] <= threshold) (*curpeak)++;
00091   if (*curpeak >= nmax) return(0); /* none found */
00092   *packpb = (*curpeak)++; /* 1st peak of current pack */
00093   npeak   = 1;            /* nb of peak in current pack */
00094 
00095   /* Find the last valid peak of the pack, that is the one before
00096      the inter-pack gap (of width FMX_TOOFAR*SNIFS_INTERSP) */
00097   while (*curpeak < nmax && 
00098          line[imax[*curpeak]] > threshold &&
00099          imax[*curpeak]-imax[(*curpeak)-1] < FMX_TOOFAR*SNIFS_INTERSP) {
00100     (*curpeak)++;
00101     npeak ++;
00102   }
00103   *packpe = (*curpeak)-1; /* last peak of current pack */
00104 
00105   return(npeak);
00106 }
00107 
00108 /* === Doxygen Comment ======================================= */
00122 /* =========================================================== */
00123 
00124 int locate_peak(max_lines line, double x)
00125 {
00126 #define x(i) (line.maxima[i].xcoord)
00127 
00128   int im;
00129   int il = -1;
00130   int iu = line.nb_max;
00131   int order = (x(iu-1) >= x(0));  /* Ascending or descending order? */
00132 
00133   while (iu-il > 1) {
00134     im = (il+iu)/2;
00135     if ((x >= x(im)) == order) il = im;
00136     else iu = im;
00137   }
00138   
00139   return(il);
00140 
00141 #undef x
00142 }
00143 
00144 /* === Doxygen Comment ======================================= */
00156 /* =========================================================== */
00157 
00158 int find_nearest_peak(Maxima_Set *maxset, int iline, float x)
00159 {
00160   int im;
00161   float xm, xp;
00162   
00163   /* Locate i such that maxima[i].xcoord <= x < maxima[i+1].xcoord */
00164   im = locate_peak(maxset->line[iline],x);
00165 
00166   /* Look for nearest peak */
00167   if (im == -1) im = 0;                      /* On the left of 1st peak */
00168   else if (im < maxset->line[iline].nb_max-1) {
00169     xm = maxset->line[iline].maxima[im  ].xcoord;
00170     xp = maxset->line[iline].maxima[im+1].xcoord;
00171     if (ABS(x - xp) < ABS(x - xm)) im++;
00172   }
00173 
00174   return(im);
00175 }
00176 
00177 /* === Doxygen Comment ======================================= */
00188 /* =========================================================== */
00189 
00190 int find_nearest_max(Maxima_Set *maxset, float x,float y, 
00191                      int *iline,int *ipeak)
00192 {
00193   float ymin,ymax;
00194   int dline;
00195 
00196   ymin = maxset->line[0].ycoord;
00197   ymax = maxset->line[maxset->nb_ycoords-1].ycoord;
00198   dline= (int)(ymax - ymin)/(maxset->nb_ycoords-1);
00199 
00200   *iline = rint((y - ymin)/dline); /* find nearest max line */
00201 
00202   if (*iline < 0 || *iline >= maxset->nb_ycoords ||
00203       maxset->line[*iline].nb_max == 0) return(UNKNOWN);
00204 
00205   /* Find closest peak in current line */
00206   *ipeak = find_nearest_peak(maxset,*iline,x);
00207 
00208   return(OK);
00209 }
00210 
00211 /* === Doxygen Comment ======================================= */
00219 /* =========================================================== */
00220 
00221 int open_max(char *maxname, Maxima_Set *maxset)
00222 {
00223   print_msg("Opening max frame '%s'...", maxname);
00224 
00225   if (load_TIGER_max(maxset,maxname)) {
00226     print_error("Cannot open maxima file %s",maxname);
00227     return(ERR_OPEN);
00228   }
00229 
00230   print_msg("   %d profiles [%.1f-%.1f]", maxset->nb_ycoords,
00231             maxset->line[0].ycoord,maxset->line[maxset->nb_ycoords-1].ycoord);
00232 
00233   return(OK);
00234 }
00235 
00236 /* === Doxygen Comment ======================================= */
00259 /* =========================================================== */
00260 
00261 int pup_get_maxdata(Maxima_Set *maxset, float xpup, float ypup, float lbdapup, 
00262                     SnifsOptics *optics, float pixsize, double lstep,
00263                     double *lbda[], double *dx[], double *sigma[])
00264 {
00265   int nlbda, nmax;
00266   int iline, ipeak;
00267   
00268   double xmla,ymla, xccd,yccd;
00269   double lambda;
00270   double maxdx = 0.3*SNIFS_INTERSP; /* Don't look further for a max */
00271   double maxsig = (optics->channel==BLUE_CHANNEL)? MAX_XDISPSIG_B:MAX_XDISPSIG_R;
00272   double minsig = (optics->channel==BLUE_CHANNEL)? MIN_XDISPSIG_B:MIN_XDISPSIG_R;
00273   
00274   /* [CCD/px] -> [MLA/mm] */
00275   snifs_optics_CCD2MLA(optics, xpup,ypup, &xmla,&ymla,pixsize,lbdapup);
00276 
00277   /* Memory allocations (maximum size) */
00278   nlbda = (optics->filter.sup_util - optics->filter.inf_util)/lstep + 1;
00279   *lbda = (double *)malloc(nlbda * sizeof(double));
00280   *dx   = (double *)malloc(nlbda * sizeof(double));
00281   *sigma= (double *)malloc(nlbda * sizeof(double));
00282 
00283   nmax = 0;
00284   for (lambda  = optics->filter.inf_util; 
00285        lambda <= optics->filter.sup_util; 
00286        lambda += lstep) { /* loop over lambda */
00287 
00288     /* [MLA/mm] -> [CCD/px] */
00289     snifs_optics_MLA2CCD(optics, xmla,ymla, &xccd,&yccd,pixsize,lambda,1); 
00290 
00291     if (!find_nearest_max(maxset, xccd,yccd, &iline,&ipeak) &&
00292         ABS(maxset->line[iline].maxima[ipeak].xcoord - xccd) < maxdx) {
00293       /* There's a max really close to there */
00294       if (maxset->line[iline].maxima[ipeak].sigma[0] > minsig &&
00295           maxset->line[iline].maxima[ipeak].sigma[0] < maxsig) {
00296         (*lbda)[nmax] = lambda;
00297         (*dx)[nmax]   = maxset->line[iline].maxima[ipeak].xcoord - xccd;
00298         (*sigma)[nmax]= maxset->line[iline].maxima[ipeak].sigma[0];
00299         nmax ++;
00300       }
00301       else if (DEBUG) {
00302         print_warning("Excluding l=%.2f: X-disp. sig = %f [%.2f,%.2f]",
00303                       lambda,maxset->line[iline].maxima[ipeak].sigma[0],
00304                       minsig,maxsig);
00305       }
00306     }
00307   }
00308 
00309   /* Reallocation (minimal size or free) */
00310   *lbda  = (double *)realloc(*lbda, nmax*sizeof(double));
00311   *dx    = (double *)realloc(*dx,   nmax*sizeof(double));
00312   *sigma = (double *)realloc(*sigma,nmax*sizeof(double));
00313 
00314   return(nmax);
00315 }
00316 

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