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 ======================================= */
00034 /* =========================================================== */
00035 
00036 long nllsqfit_2DgaussB(long *mode,long *npts,long *npar,long *ldfj,double *par,
00037                        double *f, double *fjac, long *nstate, long *iuser,
00038                        double *user)
00039 {
00040   int nx,ny, i,j,k;
00041   double dx, dy, expo, Iexp, Sx2, Sy2, dx2, dy2;
00042 
00043   nx = (int)user[0]; /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00044   ny = (int)user[1];
00045 
00046   Sx2 = SQ(par[3]);
00047   Sy2 = SQ(par[4]);
00048 
00049   for (k=0, j=0; j<ny; j++) {                /* Natural order for images */
00050     for (i=0; i<nx; i++, k++) {
00051       dx = i - par[0];
00052       dy = j - par[1];
00053       dx2 = SQ(dx);
00054       dy2 = SQ(dy);
00055       expo = exp(-(dx2/Sx2 + dy2/Sy2)/2);
00056       Iexp = par[2]*expo;
00057       f[k] = Iexp;
00058 
00059       if (*mode) {
00060         /* Partial derivatives with respect to... */
00061         fjac[k          ] = dx * f[k] / Sx2;       /* x */
00062         fjac[k + *npts  ] = dy * f[k] / Sy2;       /* y */
00063         fjac[k + *npts*2] = expo;                  /* I */
00064         fjac[k + *npts*3] = dx2*f[k]/(par[3]*Sx2); /* Sx */
00065         fjac[k + *npts*4] = dy2*f[k]/(par[4]*Sy2); /* Sy */
00066         fjac[k + *npts*5] = 1;                     /* bkgnd */
00067       }
00068       
00069       f[k] += par[5]; /* Add background */
00070     }
00071   }
00072 
00073   return(OK);
00074 }
00075 
00076 /* === Doxygen Comment ======================================= */
00093 /* =========================================================== */
00094 
00095 int main(int argc, char **argv)
00096 {
00097   int const npar = 6;        /* Parameters: (x,y) + I + (Sx,Sy) + bkgnd = 6 */
00098 
00099   IMAGE2D inframe, dbgframe;
00100 
00101   char **argval, **arglabel;
00102   char *inname, *devname;
00103 
00104   int npix[2];
00105   int i,j,k, npts, kmin, kmax, imax, jmax, status;
00106   int NOBKGND, ASCII;
00107 
00108   float min, max, xtmp, ytmp;
00109 
00110   double *x, *y, *yf;
00111   double par[npar], epar[npar], parlb[npar], parub[npar];
00112   double start[2], step[2];
00113   double rms, gof;
00114 
00115 #ifdef HAVE_LIBDISLIN
00116   Plot plot;
00117   char title[4*lg_ident+1];
00118   int nx,ny;
00119   float mincut,maxcut,**zmat;
00120 #endif
00121 
00122   set_purpose("Center a 2D-gaussian on input-frame");
00123   set_arglist("-in none -nobkgnd -dev ASCII -cuts 15,99");
00124 
00125   init_snifs("$Name:  $");
00126   init_session(argv,argc,&arglabel,&argval);
00127 
00128   if (DEBUG) {
00129     print_msg("$Id: center_gauss.c,v 1.5 2004/11/09 18:01:26 ycopin Exp $");
00130     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00131   }
00132 
00133   inname = argval[0];                        /* Input file */
00134   NOBKGND = is_true(argval[1]);              /* No background */
00135 
00136   devname = argval[2];                       /* Output device */
00137   upper_strg(devname);
00138   ASCII = !strcmp(devname,"ASCII");
00139 
00140 #ifdef HAVE_LIBDISLIN
00141   plot.Graph = !ASCII;                      /* Store flag in Plot structure */
00142   strncpy(plot.devname,devname,5);          /* Output device */
00143   if (sscanf(argval[3], "%f,%f",&mincut,&maxcut) != 2) { 
00144     print_error("Invalid option -cuts min,max (has '%s')",
00145                 argval[3]);                 /* Cuts [percentiles] */
00146     exit_session(ERR_BAD_PARAM);
00147   }
00148 #else
00149   if (!ASCII) {
00150     print_error("DISLIN output not supported in current executable");
00151     exit_session(ERR_BAD_PARAM);
00152   }
00153 #endif
00154 
00155   /* ===== ACTION! ============================== */
00156 
00157   /* ----- Open image and check type ------------------------------ */
00158 
00159   if (open_image(&inframe,inname,"I",0,"unknown"))
00160     exit_session(ERR_OPEN);
00161 
00162   if (inframe.data_type != FLOAT) 
00163     exit_session(ERR_BAD_TYPE);
00164 
00165 #ifdef HAVE_LIBDISLIN
00166 
00167   /* ----- Graphical output ------------------------------ */
00168 
00169   if (!ASCII) {
00170     plot.Limits.xmin = inframe.startx;
00171     plot.Limits.xmax = inframe.endx;
00172     plot.Limits.ymin = inframe.starty;
00173     plot.Limits.ymax = inframe.endy;
00174 
00175     plot_set_bkgnd(&inframe, &plot, mincut,maxcut, &nx,&ny, &zmat);
00176     plot_start(&plot,-0.8);                  /* Enforce isoscale */
00177     autres(nx,ny);
00178     nobar();                                 /* no color bar */
00179     plot_axes(&plot,"x","y");
00180     sursze(plot.Limits.xmin,plot.Limits.xmax,plot.Limits.ymin,plot.Limits.ymax);
00181 
00182     crvmat(zmat[0],nx,ny,1,1);
00183     //    cross();                                 /* Add central cross */
00184     sprintf(title,"%s\nCuts: [%7.1e,%7.1e]",inname,plot.mincut,plot.maxcut);
00185 
00186     free2d(&zmat,FLOAT);
00187   }
00188 #endif
00189 
00190   /* ----- minmax ------------------------------ */
00191 
00192   npts = inframe.nx * inframe.ny;            /* Total nb of points */
00193 
00194   minmax_f(inframe.data.f_data, npts, &min, &max, &kmin, &kmax);
00195 
00196   if (max == min) {
00197     print_error("Input image is uniformely equal to %f",min);
00198     close_frame(&inframe);
00199     exit_session(ERR_BAD_IMA);
00200   }
00201 
00202   /* Find position of max with modulos, since
00203      RD_frame(frame,i,j) = (frame)->data.s_data[(i)+(j)*(frame)->nx] */
00204   imax = kmax % inframe.nx;
00205   jmax = kmax / inframe.nx;
00206 
00207   if (DEBUG) {
00208     int imin,jmin;
00209     
00210     coord_frame(&inframe,imax,jmax,&xtmp,&ytmp);
00211     print_msg("   Max = %f at (%.2f x %.2f) [%d x %d]",
00212               max,xtmp,ytmp,imax,jmax);
00213 
00214     imin = kmin % inframe.nx;
00215     jmin = kmin / inframe.nx;
00216     coord_frame(&inframe,imin,jmin,&xtmp,&ytmp);
00217     print_msg("   Min = %f at (%.2f x %.2f) [%d x %d]",
00218               min,xtmp,ytmp,imin,jmin);
00219   }
00220 
00221   /* ----- 2D-gaussian fit ------------------------------ */
00222 
00223   x = (double *)calloc(npts,sizeof(double)); /* Initialize everything for dirty trick */
00224   y = (double *)malloc(npts*sizeof(double));
00225   yf= (double *)malloc(npts*sizeof(double));
00226 
00227   x[0] = inframe.nx;   /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00228   x[1] = inframe.ny; 
00229 
00230   if (NOBKGND) min = 0;                      /* No background */
00231 
00232   /* Read image data and normalize to max */
00233   /* One could also do directly:
00234      for (k=0; k<npts; k++) y[k] = (inframe.data.f_data[k] - min) / (max - min);
00235      but loosing the RD_frame interface between pgm and data */
00236   for (k=0, j=0; j<inframe.ny; j++)          /* Natural order for images */
00237     for (i=0; i<inframe.nx; i++, k++) 
00238       y[k] = (RD_frame(&inframe,i,j) - min) / (max - min);
00239 
00240   close_frame(&inframe);
00241 
00242   /* ----- Initial estimates ------------------------------ */
00243 
00244   par[0] = imax;                             /* position [px] */
00245   par[1] = jmax;
00246   par[2] = 1;                                /* normalized intensity */
00247   par[3] = 1;                                /* dispersion [px] */
00248   par[4] = 1;
00249   par[5] = 0;                                /* normalized background */
00250 
00251   /* ----- Limits ------------------------------ */
00252 
00253   parlb[0] = 0;         parub[0] = inframe.nx;    /* x,y [px] */
00254   parlb[1] = 0;         parub[1] = inframe.ny;    /* Limit to the fitting window */
00255   parlb[2] = 1e-6;      parub[2] = MAXFLOAT;      /* I > 0 */
00256   parlb[3] = 1e-6;      parub[3] = inframe.nx;
00257   parlb[4] = 1e-6;      parub[4] = inframe.ny;    /* Sx,Sy > 0 [px] */
00258   if (NOBKGND)
00259     parlb[5] = parub[5] = 0;                      /* No background */
00260   else {
00261     parlb[5] = -MAXFLOAT; parub[5] = MAXFLOAT;    /* Bkgnd */
00262   }
00263   
00264   if (DEBUG) {
00265     print_msg("Initial guess:");
00266     coord_frame(&inframe,(int)par[0],(int)par[1],&xtmp,&ytmp);
00267     print_msg("   Center:     %.2f x %.2f [%d x %d]",
00268               xtmp,ytmp,(int)par[0],(int)par[1]);
00269     print_msg("   Intensity:  %.2f [%.2f]",par[2]*(max - min),par[2]);
00270     print_msg("   Dispersion: %.2f x %.2f [%.2f x %.2f]",
00271               par[3]*inframe.stepx,par[4]*inframe.stepy, par[3],par[4]);
00272     if (!NOBKGND) 
00273       print_msg("   Background: %.2f [%.2f]",par[5]*(max - min)+min,par[5]);
00274     else 
00275       print_warning("No constant background added");
00276   }
00277 
00278   /* ----- Fit ------------------------------ */
00279 
00280   epar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00281 
00282   status = nllsqfit_bnd(x, y, NULL, npts, par, parlb, parub, npar,
00283                         *nllsqfit_2DgaussB, yf, NULL, epar, &rms,&gof);
00284 
00285   /* ----- Conclusion ------------------------------ */
00286 
00287   rms = sqrt(rms/npts);
00288 
00289   /* Unnormalize parameters and errors */
00290   //  coord_frame(inframe, par[0], par[1], &(par[0]), &(par[1])); /* Convert to WC */
00291   /* Cannot use coord_frame with non-integer pixel coords */
00292   par[0] = inframe.startx + par[0]*inframe.stepx; /* x,y [WC] */
00293   par[1] = inframe.starty + par[1]*inframe.stepy;
00294   par[2] = par[2]*(max - min);               /* I */
00295   par[3]*= inframe.stepx;                    /* Sx,Sy [WC] */
00296   par[4]*= inframe.stepx;
00297   par[5] = par[5]*(max - min) + min;         /* Bkgnd */
00298 
00299   epar[0] *= inframe.stepx;                  /* dx,dy [WC] */
00300   epar[1] *= inframe.stepy;
00301   epar[2] *= (max - min);                    /* dI */
00302   epar[3] *= inframe.stepx;                  /* dSx,dSy [WC] */
00303   epar[4] *= inframe.stepx;
00304   epar[5] *= (max - min);                    /* dBkgnd */
00305 
00306   if (VERBOSE) {
00307     print_msg("2D-gaussian fit: status=%d, RMS=%g",status,rms);
00308     print_msg("   Center:     %.2f x %.2f (+- %.2f x %.2f)",
00309               par[0],par[1],epar[0],epar[1]);
00310     print_msg("   Intensity:  %.2f        (+- %.2f)",
00311               par[2],epar[2]);
00312     print_msg("   Dispersion: %.2f x %.2f (+- %.2f x %.2f)",
00313               par[3],par[4],epar[3],epar[4]);
00314     print_msg("   FWHM:       %.2f x %.2f (+- %.2f x %.2f)",
00315               par[3]*M_FWHM,par[4]*M_FWHM,
00316               epar[3]*M_FWHM,epar[4]*M_FWHM);
00317     if (!NOBKGND) 
00318       print_msg("   Background: %.2f        (+- %.2f)",
00319                 par[5],epar[5]);
00320     else 
00321       print_warning("No constant background added");
00322   }
00323   else {
00324     printf("%f %f %f %f\n",par[0],par[1],epar[0],epar[1]);
00325     printf("%f %f\n",par[2],epar[2]);
00326     printf("%f %f %f %f\n",par[3],par[4],epar[3],epar[4]);
00327     printf("%f %f\n",par[5],epar[5]);
00328     fflush(stdout);
00329   }
00330   
00331   if (DEBUG) {
00332     npix[0] = inframe.nx;     npix[1] = inframe.ny;
00333     start[0]= inframe.startx; start[1]= inframe.starty;
00334     step[0] = inframe.stepx;  step[1] = inframe.stepy;
00335 
00336     disable_erase_flag();
00337     print_msg("DEBUG: writing fit result in dbg_center_gauss");
00338     create_frame(&dbgframe, "dbg_center_gauss", npix, start, step, 
00339                  inframe.data_type, "2D-gaussian fit",inframe.cunit);
00340     for (k=0, j=0; j<inframe.ny; j++) 
00341       for (i=0; i<inframe.nx; i++, k++) 
00342         WR_frame(&dbgframe, i,j, yf[k]*(max - min) + min);
00343     close_frame(&dbgframe);
00344     restore_erase_flag();
00345   }
00346 
00347 #ifdef HAVE_LIBDISLIN
00348 
00349   /* ----- Graphical output ------------------------------ */
00350 
00351   if (!ASCII) {
00352     int nlev = 5;                           /* Nb of levels in contour plot */
00353     float levstart, levstep;
00354     
00355     sprintf(title,"%s\nCenter: %.2fx%.2f (\\pm %.2fx%.2f)\n"
00356             "FWHM: %.2fx%.2f (\\pm %.2fx%.2f)",
00357             title,par[0],par[1],epar[0],epar[1],
00358             par[3]*M_FWHM,par[4]*M_FWHM,epar[3]*M_FWHM,epar[4]*M_FWHM);
00359     plot_title(title);
00360 
00361     alloc2d(&zmat,inframe.nx,inframe.ny,FLOAT);
00362     for (k=0, j=0; j<inframe.ny; j++)        /* Beware the order */
00363       for (i=0; i<inframe.nx; i++, k++) 
00364         zmat[i][j] = yf[k]*(max - min) + min;
00365   
00366     labels("FLOAT","CONTUR");
00367     nlev = plot_compticks(nlev, min, max, &levstart, &levstep);
00368     for (k=0; k<nlev; k++) 
00369       conmat(zmat[0],inframe.nx,inframe.ny,levstart + k*levstep);
00370 
00371     disfin();
00372   }
00373   
00374 #endif
00375 
00376   free(x); free(y); free(yf);
00377 
00378   exit_session(status);
00379   return(status);
00380 }

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