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

focus/source/center_gauss.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 
00024 #ifdef HAVE_LIBDISLIN
00025 #include <dislin.h>
00026 #endif
00027 
00028 /* === Doxygen Comment ======================================= */
00032 /* =========================================================== */
00033 
00034 long nllsqfit_2Dgauss_ima(long *mode,long *npts,long *npar,long *ldfj,double *par,
00035                           double *f, double *fjac, long *nstate, long *iuser,
00036                           double *user)
00037 {
00038   int nx,ny, i,j,k;
00039   double dx, dy, expo, Iexp, Sx2, Sy2, dx2, dy2;
00040 
00041   nx = (int)user[0]; /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00042   ny = (int)user[1];
00043 
00044   Sx2 = SQ(par[3]);
00045   Sy2 = SQ(par[4]);
00046 
00047   for (k=0, j=0; j<ny; j++) {                /* Natural order for images */
00048     for (i=0; i<nx; i++, k++) {
00049       dx = i - par[0];
00050       dy = j - par[1];
00051       dx2 = SQ(dx);
00052       dy2 = SQ(dy);
00053       expo = exp(-(dx2/Sx2 + dy2/Sy2)/2);
00054       Iexp = par[2]*expo;
00055       f[k] = Iexp;
00056 
00057       if (*mode) {
00058         /* Partial derivatives with respect to... */
00059         fjac[k          ] = dx * f[k] / Sx2;       /* x */
00060         fjac[k + *npts  ] = dy * f[k] / Sy2;       /* y */
00061         fjac[k + *npts*2] = expo;                  /* I */
00062         fjac[k + *npts*3] = dx2*f[k]/(par[3]*Sx2); /* Sx */
00063         fjac[k + *npts*4] = dy2*f[k]/(par[4]*Sy2); /* Sy */
00064         fjac[k + *npts*5] = 1;                     /* bkgnd */
00065       }
00066       
00067       f[k] += par[5]; /* Add background */
00068     }
00069   }
00070 
00071   return(OK);
00072 }
00073 
00074 /* === Doxygen Comment ======================================= */
00078 /* =========================================================== */
00079 
00080 long nllsqfit_2Daxigauss_ima(long *mode,long *npts,long *npar,long *ldfj,double *par,
00081                              double *f, double *fjac, long *nstate, long *iuser,
00082                              double *user)
00083 {
00084   int nx,ny, i,j,k;
00085   double dx, dy, expo, Iexp, S2, dx2, dy2;
00086 
00087   nx = (int)user[0]; /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00088   ny = (int)user[1];
00089 
00090   S2 = SQ(par[3]);
00091 
00092   for (k=0, j=0; j<ny; j++) {                /* Natural order for images */
00093     for (i=0; i<nx; i++, k++) {
00094       dx = i - par[0];
00095       dy = j - par[1];
00096       dx2 = SQ(dx);
00097       dy2 = SQ(dy);
00098       expo = exp(-(dx2 + dy2)/2/S2);
00099       Iexp = par[2]*expo;
00100       f[k] = Iexp;
00101 
00102       if (*mode) {
00103         /* Partial derivatives with respect to... */
00104         fjac[k          ] = dx * f[k] / S2;             /* x */
00105         fjac[k + *npts  ] = dy * f[k] / S2;             /* y */
00106         fjac[k + *npts*2] = expo;                       /* I */
00107         fjac[k + *npts*3] = (dx2+dy2)*f[k]/(par[3]*S2); /* S */
00108         fjac[k + *npts*4] = 1;                          /* bkgnd */
00109       }
00110       
00111       f[k] += par[4]; /* Add background */
00112     }
00113   }
00114 
00115   return(OK);
00116 }
00117 
00118 /* === Doxygen Comment ======================================= */
00158 /* =========================================================== */
00159 
00160 int main(int argc, char **argv)
00161 {
00162   IMAGE2D inframe, dbgframe;
00163 
00164   char **argval, **arglabel;
00165   char *inname, *devname;
00166 
00167   int npix[2];
00168   int i,j,k, npts, kmin, kmax, imax, jmax, status;
00169   int hsize, iinf, jinf, isup, jsup, inpts, jnpts;
00170   int NOBKGND, ASCII, CENTER, NOFIT, CERCLE, AXISYM;
00171   long npar;
00172   
00173   float min, max, xtmp, ytmp, xcen,ycen, xcirc,ycirc,rcirc;
00174 
00175   double *x, *y, *yf;
00176   double *par, *epar, *parlb, *parub;
00177   double start[2], step[2];
00178   double rms, gof;
00179 
00180   long (*fitting_fn)(long *mode,long *npts,long *npar,long *ldfj,double *par,
00181                      double *f, double *fjac, long *nstate, long *iuser,
00182                      double *user);
00183 
00184 #ifdef HAVE_LIBDISLIN
00185   Plot plot;
00186   char title[4*lg_ident+1];
00187   int nx,ny, nr, nw, nh;
00188   float mincut,maxcut,**zmat;
00189 #endif
00190 
00191   set_purpose("Center a 2D-gaussian on input-frame");
00192   set_arglist("-in none -nobkgnd -dev ASCII -cuts 15,99 -center NULL "
00193               "-nofit -circle NULL -axisym -hsize 0");
00194 
00195   init_snifs("$Name:  $","$Revision: 1.18 $");
00196   init_session(argv,argc,&arglabel,&argval);
00197 
00198   if (DEBUG) {
00199     print_msg("$Id: center_gauss.c,v 1.18 2005/09/15 21:43:14 ycopin Exp $");
00200     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00201   }
00202 
00203   inname = argval[0];                        /* Input file */
00204   NOBKGND = is_true(argval[1]);              /* No background */
00205 
00206   if (!strcasecmp((devname = argval[2]),"ASCII") ) /* Output device or file */
00207     *devname = '\0';                         /* Empty string = ASCII */
00208   ASCII = !strlen(devname);
00209 #ifdef HAVE_LIBDISLIN
00210   strcpy(plot.devname,devname);
00211   if (sscanf(argval[3], "%f,%f",&mincut,&maxcut) != 2) { 
00212     print_error("Invalid option -cuts min,max (has '%s')",
00213                 argval[3]);                 /* Cuts [percentiles] */
00214     exit_session(ERR_BAD_PARAM);
00215   }
00216 #else
00217   if (!ASCII) {
00218     print_error("DISLIN output not supported in current executable");
00219     exit_session(ERR_BAD_PARAM);
00220   }
00221 #endif
00222 
00223   if ((CENTER = is_set(argval[4]))) {
00224     if (sscanf(argval[4], "%f,%f",&xcen,&ycen) != 2) { 
00225       print_error("Invalid option -center x0,y0 (has '%s')",
00226                   argval[4]);                /* Center estimate */
00227       exit_session(ERR_BAD_PARAM);
00228     }
00229   }
00230 
00231   if ((NOFIT = is_true(argval[5])) && ASCII) { /* No fit done */
00232     print_error("Nothing to display (-nofit) on ASCII device (-dev ASCII)");
00233     exit_session(ERR_BAD_PARAM);
00234   }
00235 
00236   if ((CERCLE = is_set(argval[6]))) {
00237     if (sscanf(argval[6], "%f,%f,%f",&xcirc,&ycirc,&rcirc) != 3) { 
00238       print_error("Invalid option -circle X,Y,R (has '%s')",
00239                   argval[6]);                /* Plot circle */
00240       exit_session(ERR_BAD_PARAM);
00241     }
00242   }
00243 
00244   if (CERCLE && ASCII) {
00245     print_warning("No circle to display (-circle) on ASCII device (-dev ASCII)");
00246   }
00247 
00248   AXISYM = is_true(argval[7]);                /* Axisymmetric 2D-gaussian */
00249   get_argval(8,"%d",&hsize);                  /* Half-size of the fitting window */
00250 
00251   /* ===== ACTION! ============================== */
00252 
00253   /* ----- Open image and check type ------------------------------ */
00254 
00255   /* IO is needed for keyword storage */
00256   if (open_image(&inframe,inname,"IO",0,"unknown"))
00257     exit_session(ERR_OPEN);
00258 
00259   if (inframe.data_type != FLOAT) {
00260     print_error("Input frame %s is not a FLOAT image",inframe.name);
00261     exit_session(ERR_BAD_TYPE);
00262   }
00263 
00264   if (inframe.stepx != inframe.stepy && AXISYM) {
00265     print_error("Input frame %s has different steps along both axes "
00266                 "(x=%f, y=%f), cannot fit axisymmetric gaussian",
00267                 inframe.name, inframe.stepx, inframe.stepy);
00268     exit_session(ERR_BAD_PARAM);
00269   }
00270 
00271 #ifdef HAVE_LIBDISLIN
00272 
00273   /* ----- Graphical output ------------------------------ */
00274 
00275   if (!ASCII) {
00276     char xlab[lg_label+1], ylab[lg_label+1];
00277 
00278     /* Read axes labels from image */
00279     read_ima_axes(&inframe, xlab, ylab);
00280     print_msg("Axes labels: X='%s', Y='%s'",xlab,ylab);
00281 
00282     plot.Limits.xmin = inframe.startx;
00283     plot.Limits.xmax = inframe.endx;
00284     plot.Limits.ymin = inframe.starty;
00285     plot.Limits.ymax = inframe.endy;
00286 
00287     plot_set_bkgnd(&inframe, &plot, mincut,maxcut, &nx,&ny, &zmat);
00288     plot_start(&plot,-0.8);                  /* Enforce isoscale */
00289     autres(nx,ny);
00290     nobar();                                 /* no color bar */
00291 
00292     /* Tedious work around to get the full pixel display on the borders... */
00293     plot.Limits.xmin -= 0.5*inframe.stepx;
00294     plot.Limits.xmax += 0.5*inframe.stepx;
00295     plot.Limits.ymin -= 0.5*inframe.stepy;
00296     plot.Limits.ymax += 0.5*inframe.stepy;
00297     plot_axes(&plot,xlab,ylab);
00298     plot.Limits.xmin += 0.5*inframe.stepx;
00299     plot.Limits.xmax -= 0.5*inframe.stepx;
00300     plot.Limits.ymin += 0.5*inframe.stepy;
00301     plot.Limits.ymax -= 0.5*inframe.stepy;
00302     sursze(plot.Limits.xmin,plot.Limits.xmax,plot.Limits.ymin,plot.Limits.ymax);
00303 
00304     crvmat(zmat[0],nx,ny,1,1);
00305     //    cross();                                 /* Add central cross */
00306     sprintf(title,"%s\nCuts: [%7.1e,%7.1e]",inname,plot.mincut,plot.maxcut);
00307 
00308     free2d(&zmat,FLOAT);
00309 
00310     if (CERCLE) {
00311       /* Draw 1,2,3-sigma circles (beware of plot coordinates) */
00312       color("GREEN");
00313       rlsymb(3,xcirc,ycirc);
00314       nx = nxposn(xcirc);               /* Convert user coord to plot coord */
00315       ny = nyposn(ycirc);
00316       nr = nxposn(xcirc + rcirc) - nx;
00317       circle(nx,ny,nr);
00318       dash();
00319       nr = nxposn(xcirc + 2*rcirc) - nx;
00320       circle(nx,ny,nr);
00321       dot();
00322       nr = nxposn(xcirc + 3*rcirc) - nx;
00323       circle(nx,ny,nr);
00324       color("FORE");
00325       solid();
00326     }
00327 
00328     if (NOFIT) {                  /* Nothing else to be done: conclude now! */
00329       print_msg("No fit required, concluding now");
00330       plot_title(title);
00331       disfin();
00332 
00333       exit_session(OK);
00334       return(OK);
00335     }
00336 
00337     if (CENTER) {                     /* Add a red cross to center estimate */
00338       color("RED");
00339       rlsymb(3,xcen,ycen);
00340       color("FORE");
00341     }
00342   }
00343 #endif
00344 
00345   /* ----- minmax ------------------------------ */
00346 
00347   minmax_f(inframe.data.f_data, inframe.nx*inframe.ny, &min, &max, &kmin, &kmax);
00348 
00349   if (max == min) {
00350     print_error("Input image is uniformely equal to %f",min);
00351     close_frame(&inframe);
00352     exit_session(ERR_BAD_IMA);
00353   }
00354 
00355   /* ----- centroid estimate ------------------------------ */
00356 
00357   if (!CENTER) {
00358     /* Find position of max with modulos, since
00359        RD_frame(frame,i,j) = (frame)->data.s_data[(i)+(j)*(frame)->nx] */
00360     imax = kmax % inframe.nx;
00361     jmax = kmax / inframe.nx;
00362 
00363     /* A better estimate of the centroid would be to compute the centroid of
00364        the image. */
00365   }
00366   else {
00367     pixel_frame(&inframe, xcen, ycen, &imax, &jmax);
00368     max = RD_frame(&inframe,imax,jmax);
00369 
00370     print_msg("Centroid estimate: %.2f x %.2f [%d x %d]",
00371               xcen,ycen,imax,jmax);
00372   }
00373 
00374   if (DEBUG) {
00375     int imin,jmin;
00376     
00377     coord_frame(&inframe,imax,jmax,&xtmp,&ytmp);
00378     print_msg("   Max = %f at (%.2f x %.2f) [%d x %d]",
00379               max,xtmp,ytmp,imax,jmax);
00380 
00381     imin = kmin % inframe.nx;
00382     jmin = kmin / inframe.nx;
00383     coord_frame(&inframe,imin,jmin,&xtmp,&ytmp);
00384     print_msg("   Min = %f at (%.2f x %.2f) [%d x %d]",
00385               min,xtmp,ytmp,imin,jmin);
00386   }
00387 
00388   /* ----- Fitting window ------------------------------ */
00389 
00390   if (hsize) {
00391     iinf = MAX(imax - hsize, 0);
00392     jinf = MAX(jmax - hsize, 0);
00393     isup = MIN(imax + hsize, inframe.nx-1);
00394     jsup = MIN(jmax + hsize, inframe.ny-1);
00395   }
00396   else {
00397     iinf = 0;
00398     jinf = 0;
00399     isup = inframe.nx-1;
00400     jsup = inframe.ny-1;
00401   }
00402   inpts = isup - iinf + 1;
00403   jnpts = jsup - jinf + 1;
00404   npts  = inpts * jnpts;                     /* Total nb of points */
00405 
00406   /* ----- 2D-gaussian fit ------------------------------ */
00407 
00408   x = (double *)calloc(npts,sizeof(double)); /* Init. everything for dirty trick */
00409   y = (double *)malloc(npts*sizeof(double));
00410   yf= (double *)malloc(npts*sizeof(double));
00411 
00412   x[0] = inpts; /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00413   x[1] = jnpts; 
00414 
00415   if (NOBKGND) min = 0;                      /* No background */
00416 
00417   /* Read image data and normalize to max */
00418   /* One could also do directly:
00419      for (k=0; k<npts; k++) y[k] = (inframe.data.f_data[k] - min) / (max - min);
00420      but loosing the RD_frame interface between pgm and data */
00421   for (k=0, j=jinf; j<=jsup; j++)            /* Natural order for images */
00422     for (i=iinf; i<=isup; i++, k++) 
00423       y[k] = (RD_frame(&inframe,i,j) - min) / (max - min);
00424 
00425   if (AXISYM) {
00426     npar = 5;
00427     fitting_fn = nllsqfit_2Daxigauss_ima;
00428   }
00429   else {
00430     npar = 6;
00431     fitting_fn = nllsqfit_2Dgauss_ima;
00432   }
00433 
00434   par   = (double *)malloc(npar*sizeof(double));
00435   epar  = (double *)malloc(npar*sizeof(double));
00436   parlb = (double *)malloc(npar*sizeof(double));
00437   parub = (double *)malloc(npar*sizeof(double));
00438 
00439   /* ----- Initial estimates ------------------------------ */
00440 
00441   par[0] = imax - iinf;                      /* position [px] */
00442   par[1] = jmax - jinf;
00443   par[2] = 1;                                /* normalized intensity */
00444   par[3] = 1;                                /* dispersion [px] */
00445   if (!AXISYM) par[4] = 1;
00446   par[AXISYM? 4:5] = 0;                      /* normalized background */
00447 
00448   /* ----- Limits ------------------------------ */
00449 
00450   parlb[0] = 0;       parub[0] = inpts;         /* x,y [px] */
00451   parlb[1] = 0;       parub[1] = jnpts;         /* Limit to the fitting window */
00452   parlb[2] = 1e-6;    parub[2] = MAXFLOAT;      /* I > 0 */
00453   parlb[3] = 1e-6;    parub[3] = inpts;         /* dispersion > 0 [px] */
00454   if (!AXISYM) { parlb[4] = 1e-6; parub[4] = jnpts; }
00455   if (NOBKGND)                                  /* No background */
00456     parlb[AXISYM? 4:5] = parub[AXISYM? 4:5] = 0;
00457   else {                                        /* Background */
00458     parlb[AXISYM? 4:5] = -MAXFLOAT; 
00459     parub[AXISYM? 4:5] = MAXFLOAT; 
00460   } 
00461   
00462   if (DEBUG) {
00463     print_msg("Initial guess:");
00464     coord_frame(&inframe,(int)par[0]+iinf,(int)par[1]+jinf,&xtmp,&ytmp);
00465     print_msg("   Center:     %.2f x %.2f [%d x %d]",
00466               xtmp,ytmp,(int)par[0]+iinf,(int)par[1]+jinf);
00467     print_msg("   Intensity:  %.2f [%.2f]",par[2]*(max - min),par[2]);
00468     if (AXISYM)
00469       print_msg("   Dispersion: %.2f [%.2f]",
00470                 par[3]*inframe.stepx, par[3]);
00471     else
00472       print_msg("   Dispersion: %.2f x %.2f [%.2f x %.2f]",
00473                 par[3]*inframe.stepx,par[4]*inframe.stepy, par[3],par[4]);
00474     if (!NOBKGND) 
00475       print_msg("   Background: %.2f [%.2f]",
00476                 par[AXISYM? 4:5]*(max - min)+min,par[AXISYM? 4:5]);
00477     else 
00478       print_warning("No constant background added");
00479   }
00480 
00481   /* ----- Fit ------------------------------ */
00482 
00483   epar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00484 
00485   status = nllsqfit_bnd(x, y, NULL, npts, par, parlb, parub, npar,
00486                         *fitting_fn, yf, NULL, epar, &rms,&gof);
00487 
00488   /* ----- Conclusion ------------------------------ */
00489 
00490   rms = sqrt(rms/npts);
00491 
00492   /* Unnormalize parameters and errors */
00493   //  coord_frame(inframe, par[0], par[1], &(par[0]), &(par[1])); /* Convert to WC */
00494   /* Cannot use coord_frame with non-integer pixel coords */
00495   par[0] = inframe.startx + (par[0]+iinf)*inframe.stepx; /* x,y [WC] */
00496   par[1] = inframe.starty + (par[1]+jinf)*inframe.stepy;
00497   par[2] = par[2]*(max - min);               /* I */
00498   par[3]*= inframe.stepx;                    /* dispersion [WC] */
00499   if (!AXISYM) par[4]*= inframe.stepx;
00500   par[AXISYM? 4:5] = par[AXISYM? 4:5]*(max - min) + min; /* Background */
00501 
00502   epar[0] *= inframe.stepx;                  /* dx,dy [WC] */
00503   epar[1] *= inframe.stepy;
00504   epar[2] *= (max - min);                    /* dI */
00505   epar[3] *= inframe.stepx;                  /* ddisp [WC] */
00506   if (!AXISYM) epar[4] *= inframe.stepx;
00507   epar[AXISYM? 4:5] *= (max - min);          /* dBkgnd */
00508 
00509   if (VERBOSE) {
00510     print_msg("2D-gaussian fit: status=%d, RMS=%g",status,rms);
00511     print_msg("   Center:     %.2f x %.2f (+- %.2f x %.2f)",
00512               par[0],par[1],epar[0],epar[1]);
00513     print_msg("   Intensity:  %.2f        (+- %.2f)",
00514               par[2],epar[2]);
00515     if (AXISYM) {
00516       print_msg("   Dispersion: %.2f        (+- %.2f)",
00517                 par[3],epar[3]);
00518       print_msg("   FWHM:       %.2f        (+- %.2f)",
00519                 par[3]*M_FWHM,epar[3]*M_FWHM);
00520     }
00521     else {
00522       print_msg("   Dispersion: %.2f x %.2f (+- %.2f x %.2f)",
00523                 par[3],par[4],epar[3],epar[4]);
00524       print_msg("   FWHM:       %.2f x %.2f (+- %.2f x %.2f)",
00525                 par[3]*M_FWHM,par[4]*M_FWHM,
00526                 epar[3]*M_FWHM,epar[4]*M_FWHM);
00527     }
00528     if (!NOBKGND) 
00529       print_msg("   Background: %.2f        (+- %.2f)",
00530                 par[AXISYM? 4:5],epar[AXISYM? 4:5]);
00531     else 
00532       print_warning("No constant background added");
00533   }
00534   else {
00535     printf("%f %f %f %f\n",par[0],par[1],epar[0],epar[1]);
00536     printf("%f %f\n",par[2],epar[2]);
00537     if (AXISYM) printf("%f %f\n",par[3],epar[3]);
00538     else        printf("%f %f %f %f\n",par[3],par[4],epar[3],epar[4]);
00539     printf("%f %f\n",par[AXISYM? 4:5],epar[AXISYM? 4:5]);
00540     fflush(stdout);
00541   }
00542   
00543   WR_desc(&inframe,CENGAUS6,DOUBLE,npar,par);   /* Results */
00544   WR_desc(&inframe,ERRGAUS6,DOUBLE,npar,epar);  /* Errors */
00545   close_frame(&inframe);
00546 
00547   if (DEBUG) {
00548     npix[0] = inframe.nx;     npix[1] = inframe.ny;
00549     start[0]= inframe.startx; start[1]= inframe.starty;
00550     step[0] = inframe.stepx;  step[1] = inframe.stepy;
00551 
00552     disable_erase_flag();
00553     print_msg("DEBUG: writing fit result in dbg_center_gauss");
00554     create_frame(&dbgframe, "dbg_center_gauss", npix, start, step, 
00555                  inframe.data_type, "2D-gaussian fit",inframe.cunit);
00556     /* Write result in fitting window only */
00557     if (hsize)
00558       for (k=0, j=0; j<inframe.ny; j++) /* Image-oriented loop (i running faster) */
00559         for (i=0; i<inframe.nx; i++, k++) 
00560           WR_frame(&dbgframe, i,j, min);
00561     for (k=0, j=jinf; j<=jsup; j++)
00562       for (i=iinf; i<=isup; i++, k++) 
00563         WR_frame(&dbgframe, i,j, yf[k]*(max - min) + min);
00564     close_frame(&dbgframe);
00565     restore_erase_flag();
00566   }
00567 
00568 #ifdef HAVE_LIBDISLIN
00569 
00570   /* ----- Graphical output ------------------------------ */
00571 
00572   if (!ASCII) {
00573     int nlev = 5;                           /* Nb of levels in contour plot */
00574     float levstart, levstep, sig;
00575     
00576     if (hsize) {
00577       /* Draw fitting window */
00578       color("RED");
00579       coord_frame(&inframe,iinf,jsup,&xtmp,&ytmp);
00580       nx = nxposn(xtmp-.5*inframe.stepx); /* Convert user coord to plot coord */
00581       ny = nyposn(ytmp+.5*inframe.stepy); /* Beware: upper-left corner */
00582       coord_frame(&inframe,isup,jinf,&xtmp,&ytmp);
00583       nw = nxposn(xtmp+.5*inframe.stepy) - nx;
00584       nh = nyposn(ytmp-.5*inframe.stepx) - ny;
00585       rectan(nx,ny,nw,nh);
00586       color("FORE");
00587       solid();
00588     }
00589 
00590     if (AXISYM)
00591       sprintf(title,"%s\nCenter: %.2fx%.2f (\\pm %.2fx%.2f)\n"
00592               "FWHM: %.2f (\\pm %.2f)",
00593               title,par[0],par[1],epar[0],epar[1],
00594               par[3]*M_FWHM,epar[3]*M_FWHM);
00595     else
00596       sprintf(title,"%s\nCenter: %.2fx%.2f (\\pm %.2fx%.2f)\n"
00597               "FWHM: %.2fx%.2f (\\pm %.2fx%.2f)",
00598               title,par[0],par[1],epar[0],epar[1],
00599               par[3]*M_FWHM,par[4]*M_FWHM,epar[3]*M_FWHM,epar[4]*M_FWHM);
00600     plot_title(title);
00601 
00602     alloc2d(&zmat,inframe.nx,inframe.ny,FLOAT);
00603     if (hsize)
00604       for (k=0, j=0; j<inframe.ny; j++)        /* Beware the order */
00605         for (i=0; i<inframe.nx; i++, k++) 
00606           zmat[i][j] = min;
00607     for (k=0, j=jinf; j<=jsup; j++)
00608       for (i=iinf; i<=isup; i++, k++) 
00609         zmat[i][j] = yf[k]*(max - min) + min;
00610   
00611     labels("FLOAT","CONTUR");
00612     labdig(2, "CONTUR");         /* Add 2-digit labels to contours */
00613     nlev = plot_compticks(nlev, min, max, &levstart, &levstep);
00614     for (k=0; k<nlev; k++) 
00615       conmat(zmat[0],inframe.nx,inframe.ny,levstart + k*levstep);
00616 
00617     /* Draw adjusted 1,2,3-sigma circles (beware of plot coordinates) */
00618     color("BLUE");
00619     rlsymb(3,par[0],par[1]);
00620     if (AXISYM) sig = par[3];
00621     else        sig = hypot(par[3],par[4])/sqrt(2);
00622     nx = nxposn(par[0]);                /* Convert user coord to plot coord */
00623     ny = nyposn(par[1]);
00624     nr = nxposn(par[0] + sig) - nx;
00625     circle(nx,ny,nr);
00626     dash();
00627     nr = nxposn(par[0] + 2*sig) - nx;
00628     circle(nx,ny,nr);
00629     dot();
00630     nr = nxposn(par[0] + 3*sig) - nx;
00631     circle(nx,ny,nr);
00632     color("FORE");
00633     solid();
00634 
00635     disfin();
00636   }
00637   
00638 #endif
00639 
00640   free(x); free(y); free(yf);
00641   free(par); free(epar); free(parlb); free(parub);
00642 
00643   exit_session(status);
00644   return(status);
00645 }

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