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

focus/source/center_2gauss.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00014 /* =========================================================== */
00015 
00016 #include <IFU_io.h>
00017 #include <IFU_math.h>
00018 #include <snifs.h>
00019 #include <focus.h>
00020 #include "../../extract/incl/extract.h" /* Probably not the right way... */
00021 
00022 #ifdef HAVE_LIBDISLIN
00023 #include <dislin.h>
00024 #endif
00025 
00026 /* === Doxygen Comment ======================================= */
00031 /* =========================================================== */
00032 
00033 long nllsqfit_double_2Daxigauss_ima(long *mode,long *npts,long *npar,long *ldfj,double *par,
00034                                     double *f, double *fjac, long *nstate, long *iuser,
00035                                     double *user)
00036 {
00037   int nx,ny, i,j,k;
00038   double dx[2], dy[2], expo[2], Iexp[2], S2, dx2[2], dy2[2];
00039 
00040   nx = (int)user[0]; /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00041   ny = (int)user[1];
00042 
00043   S2 = SQ(par[6]);
00044 
00045   for (k=0, j=0; j<ny; j++) {                /* Natural order for images */
00046     for (i=0; i<nx; i++, k++) {
00047       dx[0] = i - par[0]; dy[0] = j - par[1];
00048       dx[1] = i - par[2]; dy[1] = j - par[3];
00049       dx2[0] = SQ(dx[0]); dy2[0] = SQ(dy[0]);
00050       dx2[1] = SQ(dx[1]); dy2[1] = SQ(dy[1]);
00051       expo[0] = exp(-(dx2[0] + dy2[0])/2/S2);
00052       expo[1] = exp(-(dx2[1] + dy2[1])/2/S2);
00053       Iexp[0] = par[4]*expo[0];
00054       Iexp[1] = par[5]*expo[1];
00055       f[k] = Iexp[0] + Iexp[1];
00056 
00057       if (*mode) {
00058         /* Partial derivatives with respect to... */
00059         fjac[k          ] = dx[0]/S2 * Iexp[0]; /* x1 */
00060         fjac[k + *npts  ] = dy[0]/S2 * Iexp[0]; /* y1 */
00061         fjac[k + *npts*2] = dx[1]/S2 * Iexp[1]; /* x2 */
00062         fjac[k + *npts*3] = dy[1]/S2 * Iexp[1]; /* y2 */
00063         fjac[k + *npts*4] = expo[0];            /* I1 */
00064         fjac[k + *npts*5] = expo[1];            /* I2 */
00065         fjac[k + *npts*6] = ( (dx2[0]+dy2[0])*Iexp[0] + (dx2[1]+dy2[1])*Iexp[1] ) / 
00066           (par[6]*S2);                          /* S */
00067         fjac[k + *npts*7] = 1;                  /* bkgnd */
00068       }
00069       
00070       f[k] += par[7]; /* Add background */
00071     }
00072   }
00073 
00074   return(OK);
00075 }
00076 
00077 /* === Doxygen Comment ======================================= */
00089 /* =========================================================== */
00090 
00091 int main(int argc, char **argv)
00092 {
00093   IMAGE2D inframe, dbgframe;
00094 
00095   char **argval, **arglabel;
00096   char *inname, *devname;
00097 
00098   long const npar = 8;
00099 
00100   int npix[2];
00101   int i,j,k, npts, kmin, kmax, imax[2], jmax[2], status;
00102   int ASCII;
00103   
00104   float min, max, peak[2], xtmp, ytmp, xcen[2],ycen[2];
00105 
00106   double *x, *y, *yf;
00107   double par[npar], epar[npar], parlb[npar], parub[npar];
00108   double start[2], step[2];
00109   double rms, gof, dx, dy, sep;
00110 
00111 #ifdef HAVE_LIBDISLIN
00112   Plot plot;
00113   char title[4*lg_ident+1];
00114   int nx,ny, nr;
00115   float mincut,maxcut,**zmat;
00116 #endif
00117 
00118   set_purpose("Center a double 2D axisymmetric gaussian on input-frame");
00119   set_arglist("-in none -dev ASCII -cuts 15,99 -center NONE");
00120 
00121   init_snifs("$Name:  $","$Revision: 1.3 $");
00122   init_session(argv,argc,&arglabel,&argval);
00123 
00124   if (DEBUG) {
00125     print_msg("$Id: center_2gauss.c,v 1.3 2005/09/15 21:43:14 ycopin Exp $");
00126     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00127   }
00128 
00129   inname = argval[0];                        /* Input file */
00130 
00131   if (!strcasecmp((devname = argval[1]),"ASCII") ) /* Output device or file */
00132     *devname = '\0';                         /* Empty string = ASCII */
00133   ASCII = !strlen(devname);
00134 #ifdef HAVE_LIBDISLIN
00135   strcpy(plot.devname,devname);
00136   if (sscanf(argval[2], "%f,%f",&mincut,&maxcut) != 2) { 
00137     print_error("Invalid option -cuts min,max (has '%s')",
00138                 argval[2]);                 /* Cuts [percentiles] */
00139     exit_session(ERR_BAD_PARAM);
00140   }
00141 #else
00142   if (!ASCII) {
00143     print_error("DISLIN output not supported in current executable");
00144     exit_session(ERR_BAD_PARAM);
00145   }
00146 #endif
00147 
00148   if (sscanf(argval[3], "%f,%f,%f,%f",
00149              &(xcen[0]),&(ycen[0]),&(xcen[1]),&(ycen[1])) != 4) { 
00150     print_error("Invalid option -center x1,y1,x2,y2 (has '%s')",
00151                 argval[3]);                /* Center estimates */
00152     exit_session(ERR_BAD_PARAM);
00153   }
00154 
00155   /* ===== ACTION! ============================== */
00156 
00157   /* ----- Open image and check type ------------------------------ */
00158 
00159   /* IO is needed for keyword storage */
00160   if (open_image(&inframe,inname,"IO",0,"unknown"))
00161     exit_session(ERR_OPEN);
00162 
00163   if (inframe.data_type != FLOAT) {
00164     print_error("Input frame %s is not a FLOAT image",inframe.name);
00165     exit_session(ERR_BAD_TYPE);
00166   }
00167 
00168 #ifdef HAVE_LIBDISLIN
00169 
00170   /* ----- Graphical output ------------------------------ */
00171 
00172   if (!ASCII) {
00173     char xlab[lg_label+1], ylab[lg_label+1];
00174 
00175     /* Read axes labels from image */
00176     read_ima_axes(&inframe, xlab, ylab);
00177     print_msg("Axes labels: X='%s', Y='%s'",xlab,ylab);
00178 
00179     plot.Limits.xmin = inframe.startx;
00180     plot.Limits.xmax = inframe.endx;
00181     plot.Limits.ymin = inframe.starty;
00182     plot.Limits.ymax = inframe.endy;
00183 
00184     plot_set_bkgnd(&inframe, &plot, mincut,maxcut, &nx,&ny, &zmat);
00185     plot_start(&plot,-0.8);                  /* Enforce isoscale */
00186     autres(nx,ny);
00187     nobar();                                 /* no color bar */
00188 
00189     /* Tedious work around to get the full pixel display on the borders... */
00190     plot.Limits.xmin -= 0.5*inframe.stepx;
00191     plot.Limits.xmax += 0.5*inframe.stepx;
00192     plot.Limits.ymin -= 0.5*inframe.stepy;
00193     plot.Limits.ymax += 0.5*inframe.stepy;
00194     plot_axes(&plot,xlab,ylab);
00195     plot.Limits.xmin += 0.5*inframe.stepx;
00196     plot.Limits.xmax -= 0.5*inframe.stepx;
00197     plot.Limits.ymin += 0.5*inframe.stepy;
00198     plot.Limits.ymax -= 0.5*inframe.stepy;
00199     sursze(plot.Limits.xmin,plot.Limits.xmax,plot.Limits.ymin,plot.Limits.ymax);
00200 
00201     crvmat(zmat[0],nx,ny,1,1);
00202     //    cross();                                 /* Add central cross */
00203     sprintf(title,"%s\nCuts: [%7.1e,%7.1e]",inname,plot.mincut,plot.maxcut);
00204 
00205     free2d(&zmat,FLOAT);
00206 
00207     color("RED");                     /* Add a red cross to center estimates */
00208     rlsymb(3,xcen[0],ycen[0]);
00209     rlsymb(3,xcen[1],ycen[1]);
00210     color("FORE");
00211   }
00212 #endif
00213 
00214   /* ----- minmax ------------------------------ */
00215 
00216   npts = inframe.nx * inframe.ny;            /* Total nb of points */
00217 
00218   minmax_f(inframe.data.f_data, npts, &min, &max, &kmin, &kmax);
00219 
00220   if (max == min) {
00221     print_error("Input image is uniformely equal to %f",min);
00222     close_frame(&inframe);
00223     exit_session(ERR_BAD_IMA);
00224   }
00225 
00226   /* ----- centroid estimate ------------------------------ */
00227 
00228   pixel_frame(&inframe, xcen[0], ycen[0], &(imax[0]), &(jmax[0]));
00229   peak[0] = RD_frame(&inframe,imax[0],jmax[0]);
00230 
00231   pixel_frame(&inframe, xcen[1], ycen[1], &(imax[1]), &(jmax[1]));
00232   peak[1] = RD_frame(&inframe,imax[1],jmax[1]);
00233 
00234   print_msg("Centroid estimates: #1: %.2f x %.2f [%d x %d], #2: %.2f x %.2f [%d x %d]",
00235             xcen[0],ycen[0],imax[0],jmax[0], xcen[1],ycen[1],imax[1],jmax[1]);
00236 
00237   /* ----- 2D-gaussian fit ------------------------------ */
00238 
00239   x = (double *)calloc(npts,sizeof(double)); /* Initialize everything for dirty trick */
00240   y = (double *)malloc(npts*sizeof(double));
00241   yf= (double *)malloc(npts*sizeof(double));
00242 
00243   x[0] = inframe.nx;   /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00244   x[1] = inframe.ny; 
00245 
00246   /* Read image data and normalize to max */
00247   /* One could also do directly:
00248      for (k=0; k<npts; k++) y[k] = (inframe.data.f_data[k] - min) / (max - min);
00249      but loosing the RD_frame interface between pgm and data */
00250   for (k=0, j=0; j<inframe.ny; j++)          /* Natural order for images */
00251     for (i=0; i<inframe.nx; i++, k++) 
00252       y[k] = (RD_frame(&inframe,i,j) - min) / (max - min);
00253 
00254   /* ----- Initial estimates ------------------------------ */
00255 
00256   par[0] = imax[0]; par[1] = jmax[0];        /* positions [px] */
00257   par[2] = imax[1]; par[3] = jmax[1];
00258   par[4] = 1;       par[5] = 1;              /* normalized intensities */
00259   par[6] = 1;                                /* dispersion [px] */
00260   par[7] = 0;                                /* normalised background */
00261 
00262   /* ----- Limits ------------------------------ */
00263 
00264   parlb[0] = 0;         parub[0] = inframe.nx; /* x,y [px] */
00265   parlb[1] = 0;         parub[1] = inframe.ny; /* Limit to the input frame */
00266   parlb[2] = 0;         parub[2] = inframe.nx;
00267   parlb[3] = 0;         parub[3] = inframe.ny;
00268   parlb[4] = 1e-6;      parub[4] = MAXFLOAT;   /* I > 0 */
00269   parlb[5] = 1e-6;      parub[5] = MAXFLOAT;
00270   parlb[6] = 1e-6;      parub[6] = inframe.nx; /* dispersion > 0 [px] */
00271   parlb[7] = -MAXFLOAT; parub[7] = MAXFLOAT; 
00272   
00273   if (DEBUG) {
00274     print_msg("Initial guess:");
00275     coord_frame(&inframe,(int)par[0],(int)par[1],&xtmp,&ytmp);
00276     print_msg("   Center #1:     %.2f x %.2f [%d x %d]",
00277               xtmp,ytmp,(int)par[0],(int)par[1]);
00278     coord_frame(&inframe,(int)par[2],(int)par[3],&xtmp,&ytmp);
00279     print_msg("   Center #2:     %.2f x %.2f [%d x %d]",
00280               xtmp,ytmp,(int)par[2],(int)par[3]);
00281     print_msg("   Intensity #1:  %.2f [%.2f]",par[4]*(max - min),par[4]);
00282     print_msg("   Intensity #2:  %.2f [%.2f]",par[5]*(max - min),par[5]);
00283     print_msg("   Dispersion:    %.2f [%.2f]",par[6]*inframe.stepx, par[6]);
00284     print_msg("   Background:    %.2f [%.2f]",par[7]*(max - min)+min,par[7]);
00285   }
00286 
00287   /* ----- Fit ------------------------------ */
00288 
00289   epar[0] = 2; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00290 
00291   status = nllsqfit_bnd(x, y, NULL, npts, par, parlb, parub, npar,
00292                         *nllsqfit_double_2Daxigauss_ima, yf, NULL, epar, &rms,&gof);
00293 
00294   /* ----- Conclusion ------------------------------ */
00295 
00296   rms = sqrt(rms/npts);
00297 
00298   /* Unnormalize parameters and errors */
00299   //  coord_frame(inframe, par[0], par[1], &(par[0]), &(par[1])); /* Convert to WC */
00300   /* Cannot use coord_frame with non-integer pixel coords */
00301   par[0] = inframe.startx + par[0]*inframe.stepx; /* x,y [WC] */
00302   par[1] = inframe.starty + par[1]*inframe.stepy;
00303   par[2] = inframe.startx + par[2]*inframe.stepx;
00304   par[3] = inframe.starty + par[3]*inframe.stepy;
00305   par[4] = par[4]*(max - min);               /* I */
00306   par[5] = par[5]*(max - min);
00307   par[6]*= inframe.stepx;                    /* dispersion [WC] */
00308   par[7] = par[7]*(max - min) + min;         /* Background */
00309 
00310   epar[0] *= inframe.stepx;                  /* dx,dy [WC] */
00311   epar[1] *= inframe.stepy;
00312   epar[2] *= inframe.stepx;
00313   epar[3] *= inframe.stepy;
00314   epar[4] *= (max - min);                    /* dI */
00315   epar[5] *= (max - min);
00316   epar[6] *= inframe.stepx;                  /* ddisp [WC] */
00317   epar[7] *= (max - min);                    /* dBkgnd */
00318 
00319   print_msg("Double 2D axisymmetric gaussian fit: status=%d, RMS=%g",status,rms);
00320   print_msg("   Center #1:     %.2f x %.2f (+- %.2f x %.2f)",
00321             par[0],par[1],epar[0],epar[1]);
00322   print_msg("   Center #2:     %.2f x %.2f (+- %.2f x %.2f)",
00323             par[2],par[3],epar[2],epar[3]);
00324   dx = par[2]-par[0];
00325   dy = par[3]-par[1];
00326   sep = hypot(dx,dy);
00327   print_msg("   Separation:    %.2f        (+- %.2f)",
00328             sep,
00329             sqrt( SQ(dx)*(SQ(epar[0])+SQ(epar[2])) + 
00330                   SQ(dy)*(SQ(epar[1])+SQ(epar[3])) ) / sep);
00331   print_msg("   Position angle (1->2)/up: %.2f deg (+- %.2f)",
00332             atan2(-dx, dy)/M_PI*180.,
00333             sqrt( ( (SQ(epar[0])+SQ(epar[2]))/SQ(dy) + 
00334                     (SQ(epar[1])+SQ(epar[3]))*SQ(dx/SQ(dy)) ) / 
00335                   SQ(1+SQ(dx/dy)) )/M_PI*180.);
00336   print_msg("   Intensity #1:  %.2f        (+- %.2f)",
00337             par[4],epar[4]);
00338   print_msg("   Intensity #2:  %.2f        (+- %.2f)",
00339             par[5],epar[5]);
00340   print_msg("   Dispersion:    %.2f        (+- %.2f)",
00341             par[6],epar[6]);
00342   print_msg("   FWHM:          %.2f        (+- %.2f)",
00343             par[6]*M_FWHM,epar[6]*M_FWHM);
00344   print_msg("   Background:    %.2f        (+- %.2f)",
00345             par[7],epar[7]);
00346 
00347   close_frame(&inframe);
00348 
00349   if (DEBUG) {
00350     npix[0] = inframe.nx;     npix[1] = inframe.ny;
00351     start[0]= inframe.startx; start[1]= inframe.starty;
00352     step[0] = inframe.stepx;  step[1] = inframe.stepy;
00353 
00354     disable_erase_flag();
00355     print_msg("DEBUG: writing fit result in dbg_center_gauss");
00356     create_frame(&dbgframe, "dbg_center_2gauss", npix, start, step, 
00357                  inframe.data_type, "Double 2D axisym. gaussian fit",inframe.cunit);
00358     for (k=0, j=0; j<inframe.ny; j++) 
00359       for (i=0; i<inframe.nx; i++, k++) 
00360         WR_frame(&dbgframe, i,j, yf[k]*(max - min) + min);
00361     close_frame(&dbgframe);
00362     restore_erase_flag();
00363   }
00364 
00365 #ifdef HAVE_LIBDISLIN
00366 
00367   /* ----- Graphical output ------------------------------ */
00368 
00369   if (!ASCII) {
00370     int nlev = 5;                           /* Nb of levels in contour plot */
00371     float levstart, levstep, sig;
00372     
00373     sprintf(title,"%s\nCenter: %.2fx%.2f (\\pm %.2fx%.2f), %.2fx%.2f (\\pm %.2fx%.2f)\n"
00374             "FWHM: %.2f (\\pm %.2f)",
00375             title,par[0],par[1],epar[0],epar[1],par[2],par[3],epar[2],epar[3],
00376             par[6]*M_FWHM,epar[6]*M_FWHM);
00377     plot_title(title);
00378 
00379     alloc2d(&zmat,inframe.nx,inframe.ny,FLOAT);
00380     for (k=0, j=0; j<inframe.ny; j++)        /* Beware the order */
00381       for (i=0; i<inframe.nx; i++, k++) 
00382         zmat[i][j] = yf[k]*(max - min) + min;
00383   
00384     labels("FLOAT","CONTUR");
00385     nlev = plot_compticks(nlev, min, max, &levstart, &levstep);
00386     for (k=0; k<nlev; k++) 
00387       conmat(zmat[0],inframe.nx,inframe.ny,levstart + k*levstep);
00388 
00389     /* Draw adjusted 3-sigma circles (beware of plot coordinates) */
00390     color("BLUE");
00391     rlsymb(3,par[0],par[1]);
00392     rlsymb(3,par[2],par[3]);
00393     sig = par[6];
00394     dash();
00395     nx = nxposn(par[0]);                /* Convert user coord to plot coord */
00396     ny = nyposn(par[1]);
00397     nr = nxposn(par[0] + sig) - nx;
00398     circle(nx,ny,nr);
00399     nx = nxposn(par[2]);                /* Convert user coord to plot coord */
00400     ny = nyposn(par[3]);
00401     nr = nxposn(par[2] + sig) - nx;
00402     circle(nx,ny,nr);
00403     color("FORE");
00404     solid();
00405 
00406     disfin();
00407   }
00408   
00409 #endif
00410 
00411   free(x); free(y); free(yf);
00412 
00413   exit_session(status);
00414   return(status);
00415 }

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