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

extract/lib/proc_corr.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 
00021 #define DBG_PREFIX "dbg_corr"
00022 
00024 #define TOO_CLOSE 0.1
00025 
00026 /* === Doxygen Comment ======================================= */
00035 /* =========================================================== */
00036 
00037 int extract_arc(const IMAGE2D *frame, int npix, double **tab)
00038 {
00039   int i, j, i1, i2, j1, j2;
00040 
00041   static int index = 1; /* Used to differentiate debug output */
00042   char dbgname[lg_name+1];
00043 
00044   /* Store central part of frame */
00045   if (frame->nx < npix || frame->ny < npix) {
00046     print_error("Frame %s (%dx%d) is <%d in one direction.", 
00047                 frame->name,frame->nx,frame->ny, npix);
00048     return(ERR_BAD_SIZE);
00049   }
00050 
00051   i1 = (frame->nx - npix)/2; 
00052   i2 = (frame->nx + npix)/2;
00053   j1 = (frame->ny - npix)/2; 
00054   j2 = (frame->ny + npix)/2;
00055   for (j=j1; j<j2; j++)
00056     for (i=i1; i<i2; i++) 
00057       tab[i-i1][j-j1] = RD_frame(frame,i,j);
00058 
00059   if (DEBUG) {
00060     int npixdbg[2];
00061 
00062     npixdbg[0] = i2 - i1; 
00063     npixdbg[1] = j2 - j1;
00064 
00065     sprintf(dbgname, "%s_arc%0d", DBG_PREFIX,index);
00066     print_msg("DEBUG: Storing central part (%dx%d) of %s in %s", 
00067               npixdbg[0], npixdbg[1], frame->name, dbgname);
00068 
00069     disable_erase_flag();
00070     dump_2Darray((void **)tab,DOUBLE,npixdbg,dbgname);
00071     restore_erase_flag();
00072 
00073     index ++;
00074   }
00075 
00076   return(OK);
00077 }
00078 
00079 /* === Doxygen Comment ======================================= */
00085 /* =========================================================== */
00086 
00087 void negate_2D(double **I, int nx, int ny)
00088 {
00089   int k,l;
00090 
00091   for (k=0; k<nx; k++)
00092     for (l=0; l<ny; l++) 
00093       I[k][l] = -I[k][l];
00094 }
00095 
00096 /* === Doxygen Comment ======================================= */
00109 /* =========================================================== */
00110 
00111 void complex_product_2D(double **AR, double **AI, 
00112                         double **BR, double **BI, 
00113                         int nx, int ny, 
00114                         double **CR, double **CI)
00115 {
00116   int i,j;
00117   
00118   for (i=0; i<nx; i++)
00119     for (j=0; j<ny; j++) {
00120       CR[i][j] = AR[i][j]*BR[i][j] - AI[i][j]*BI[i][j];
00121       CI[i][j] = AR[i][j]*BI[i][j] + AI[i][j]*BR[i][j];
00122     }
00123 }
00124 
00125 /* === Doxygen Comment ======================================= */
00138 /* =========================================================== */
00139 
00140 int FFT2D_real(double **in, int nx, int ny, 
00141                double **outR, double **outI)
00142 {
00143   char init='i';
00144   int i,j;
00145   long m,n,ifail;
00146   double *trigm, *trign, *work;
00147   long c06fuf();
00148   
00149   /* Working space */
00150   trigm = (double *)malloc(2*nx   *sizeof(double));
00151   trign = (double *)malloc(2*ny   *sizeof(double));
00152   work  = (double *)malloc(2*nx*ny*sizeof(double));
00153 
00154   /* Dump input real array */
00155   for (i=0; i<nx; i++)
00156     for (j=0; j<ny; j++) {
00157       outR[i][j] = in[i][j];
00158       outI[i][j] = 0;
00159     }
00160 
00161   /* FFT (initial call) */
00162   if (DEBUG) ifail = -1; /* noisy soft fail */
00163   else       ifail =  1; /* silent soft fail */
00164   m = nx;
00165   n = ny;
00166   c06fuf(&n, &m, outR[0], outI[0], &init, trigm, trign, work, &ifail);
00167 
00168   /* Memory management */
00169   free(trigm); free(trign); free(work);
00170 
00171   if (ifail) /* Should never happen... */
00172     print_error("Status = %ld in FFT2D_real (%dx%d)",ifail,nx,ny);
00173 
00174   return((int)ifail); 
00175 }
00176 
00177 /* === Doxygen Comment ======================================= */
00190 /* =========================================================== */
00191 
00192 int iFFT2D_real(double **inR, double **inI, 
00193                 int nx, int ny, double **outR)
00194 {
00195   char init='i';
00196   int i,j,ii,jj;
00197   long m,n,m2,n2,ifail;
00198   double *trigm, *trign, *work;
00199   double norm;
00200   long c06fuf();
00201 
00202   /* Working space */
00203   trigm = (double *)malloc(2*nx   *sizeof(double));
00204   trign = (double *)malloc(2*ny   *sizeof(double));
00205   work  = (double *)malloc(2*nx*ny*sizeof(double));
00206 
00207   /* Inverse FFT */
00208   if (DEBUG) ifail = -1; /* noisy soft fail */
00209   else       ifail =  1; /* silent soft fail */
00210   m = nx;
00211   n = ny;
00212 
00213   negate_2D(inI, nx,ny);
00214   c06fuf(&n, &m, inR[0], inI[0], &init, trigm, trign, work, &ifail);
00215 
00216   /* Memory management */
00217   free(trigm); free(trign); free(work);
00218 
00219   if (ifail) /* Should never happen... */
00220     print_error("Status = %ld in FFT2D_real (%dx%d)",ifail,nx,ny);
00221 
00222   /* Swap, normalize and save into output array */
00223   m2 = m/2; 
00224   n2 = n/2;
00225   norm = sqrt(m*n);
00226   for (i=0, ii=m2; i<m2; i++, ii++) 
00227     for (j=0, jj=n2; j<n2; j++, jj++) {
00228       outR[ii][jj]= norm*inR[i][j];    
00229       outR[i][j]  = norm*inR[m-m2+i][n-n2+j]; 
00230       outR[ii][j] = norm*inR[i][n-n2+j]; 
00231       outR[i][jj] = norm*inR[m-m2+i][j];
00232     }
00233 
00234   return((int)ifail);
00235 }
00236 
00237 /* === Doxygen Comment ======================================= */
00248 /* =========================================================== */
00249 
00250 int correlation2D(double **A, double **B, 
00251                   int nx, int ny, double **X)
00252 {
00253   long ifail;
00254   double **tr, **ti, **gr, **gi, **fr, **fi;
00255 
00256   /* Memory allocation */
00257   alloc2d(&tr,nx,ny,DOUBLE); alloc2d(&ti,nx,ny,DOUBLE); 
00258   alloc2d(&gr,nx,ny,DOUBLE); alloc2d(&gi,nx,ny,DOUBLE);
00259   alloc2d(&fr,nx,ny,DOUBLE); alloc2d(&fi,nx,ny,DOUBLE);
00260 
00261   /* FFTs of ref_tab and cur_tab */
00262   if ((ifail = FFT2D_real(A, nx,ny, tr,ti))) goto clean;
00263   if ((ifail = FFT2D_real(B, nx,ny, gr,gi))) goto clean;
00264 
00265   /* FT(correlation) = FT(T) x FT(G)* */
00266   negate_2D(gi, nx,ny);
00267   complex_product_2D(tr,ti, gr,gi, nx,ny, fr,fi);
00268 
00269   /* Inverse FFT */
00270   if ((ifail = iFFT2D_real(fr,fi, nx,ny, X))) goto clean;
00271 
00272   /* Memory management */
00273  clean:
00274   free2d(&tr,DOUBLE); free2d(&ti,DOUBLE); 
00275   free2d(&gr,DOUBLE); free2d(&gi,DOUBLE);
00276   free2d(&fr,DOUBLE); free2d(&fi,DOUBLE);
00277 
00278   return((int)ifail);
00279 }
00280 
00281 /* === Doxygen Comment ======================================= */
00299 /* =========================================================== */
00300 
00301 long nllsqfit_2Dgauss(long *mode,long *npts,long *npar,long *ldfj,double *par,
00302                       double *f, double *fjac, long *nstate, long *iuser,
00303                       double *user)
00304 {
00305   int nx,ny, i,j,k;
00306   double dx, dy, expo, Iexp, Sx2, Sy2, dx2, dy2;
00307 
00308   nx = (int)user[0]; /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00309   ny = (int)user[1];
00310 
00311   Sx2 = SQ(par[3]);
00312   Sy2 = SQ(par[4]);
00313 
00314   for (k=0, i=0; i<nx; i++) {                /* Natural order for 2D-arrays */
00315     for (j=0; j<ny; j++, k++) {
00316       dx = i - par[0];
00317       dy = j - par[1];
00318       dx2 = SQ(dx);
00319       dy2 = SQ(dy);
00320       expo = exp(-(dx2/Sx2 + dy2/Sy2)/2);
00321       Iexp = par[2]*expo;
00322       f[k] = Iexp;
00323 
00324       if (*mode) {
00325         /* Partial derivatives with respect to... */
00326         fjac[k          ] = dx * f[k] / Sx2;       /* x */
00327         fjac[k + *npts  ] = dy * f[k] / Sy2;       /* y */
00328         fjac[k + *npts*2] = expo;                  /* I */
00329         fjac[k + *npts*3] = dx2*f[k]/(par[3]*Sx2); /* Sx */
00330         fjac[k + *npts*4] = dy2*f[k]/(par[4]*Sy2); /* Sy */
00331         fjac[k + *npts*5] = 1;                     /* bkgnd */
00332       }
00333       
00334       f[k] += par[5]; /* Add background */
00335     }
00336   }
00337 
00338   return(OK);
00339 }
00340 
00341 /* === Doxygen Comment ======================================= */
00359 /* =========================================================== */
00360 
00361 int Fit_Xpeak(double ** ref_tab, double ** cur_tab, 
00362               int nx, int ny, 
00363               int winx, int winy, double par[6], double epar[6])
00364 {
00365   static int index = 1;               /* Used to differentiate debug output */
00366   char dbgname[lg_name+1];
00367   int  npts, status, npix[2], i1, i2, j1, j2;
00368   int const npar = 6;        /* Parameters: (x,y) + I + (Sx,Sy) + bkgnd = 6 */
00369   long i, j, ifail, k, imax, jmax;
00370   double **corr; 
00371   double parlb[npar], parub[npar], *x, *y, *yf;
00372   double rms, gof;
00373   double Imax, Imin;
00374 
00375   /* ===== Compute X-correlation ============================== */
00376 
00377   /* Memory allocation */
00378   alloc2d(&corr,nx,ny,DOUBLE);
00379 
00380   ifail = correlation2D(ref_tab, cur_tab, nx,ny, corr);
00381 
00382   if (DEBUG) {
00383     npix[0] = nx;
00384     npix[1] = ny;
00385 
00386     sprintf(dbgname, "%s_%0d", DBG_PREFIX,index);
00387     print_msg("DEBUG: Storing arc X-corr (%dx%d) in %s",
00388               npix[0],npix[1],dbgname);
00389     disable_erase_flag();
00390     dump_2Darray((void **)corr,DOUBLE,npix,dbgname);
00391     restore_erase_flag();
00392   }
00393 
00394   /* ===== Central peak fit ============================== */
00395 
00396   /* Fitting window limits [px] */
00397   i1 = nx/2 - winx;  i2 = nx/2 + winx - 1;
00398   j1 = ny/2 - winy;  j2 = ny/2 + winy - 1;
00399 
00400   if (DEBUG) 
00401     print_msg("DEBUG: correlation peak fitting window: [%d,%d:%d,%d]",
00402               i1,j1,i2,j2);
00403 
00404   if (i1<0 || i2>=nx || j1<0 || j2>=ny) {
00405     print_error("Fitting window larger than full window (%dx%d)",nx,ny);
00406     return(ERR_BAD_SIZE);
00407   }
00408 
00409   npts = 4*winx*winy;
00410   
00411   x = (double *)calloc(npts,sizeof(double)); /* Initialize everything for dirty trick */
00412   y = (double *)malloc(npts*sizeof(double));
00413   yf= (double *)malloc(npts*sizeof(double));
00414 
00415   x[0] = 2*winx; /* Dirty trick to pass 2D-sizes of y to nllsqfit_bnd */
00416   x[1] = 2*winy; 
00417 
00418   /* ----- Estimates ------------------------------ */
00419 
00420   /* Fill in the y-array and look for peak position [px] and intensity */
00421   Imax = -MAXFLOAT;
00422   Imin =  MAXFLOAT;
00423   for (k=0, i=i1; i<=i2; i++) {
00424     for (j=j1; j<=j2; j++, k++) {
00425       y[k] = corr[i][j];
00426 
00427       if (y[k] > Imax) {                     /* Max and abs. position */
00428         Imax = y[k];
00429         imax = i;
00430         jmax = j;
00431       }
00432       if (y[k] < Imin) Imin = y[k];          /* Min */
00433     }
00434   }
00435   imax -= i1;                                /* Max relative position */
00436   jmax -= j1;
00437 
00438   /* Normalize by Imax */
00439   for (k=0, i=i1; i<=i2; i++)
00440     for (j=j1; j<=j2; j++, k++)
00441       y[k] /= Imax;
00442   Imin /= Imax;
00443 
00444   par[0] = imax; /* position [px] */
00445   par[1] = jmax;
00446   par[2] = 1;    /* normalized intensity */
00447   par[3] = 0.8;  /* dispersion [px] */
00448   par[4] = 0.8;
00449   par[5] = Imin; /* normalized background */
00450 
00451   /* ----- Limits ------------------------------ */
00452 
00453   parlb[0] = 0;         parub[0] = 2*winx + 1; /* x,y [px] */
00454   parlb[1] = 0;         parub[1] = 2*winy + 1; /* Limit to the fitting window */
00455   parlb[2] = 0;         parub[2] = MAXFLOAT;   /* I */
00456   parlb[3] = 0.1;       parub[3] = MAXFLOAT;   /* Sx,Sy [px] */
00457   parlb[4] = 0.1;       parub[4] = MAXFLOAT;
00458   parlb[5] = 0;         parub[5] = MAXFLOAT;   /* Bkgnd */
00459 
00460   /* ----- Fit ------------------------------ */
00461 
00462   epar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00463 
00464   status = nllsqfit_bnd(x, y, NULL, npts, par, parlb, parub, npar,
00465                         *nllsqfit_2Dgauss, yf, NULL, epar, &rms,&gof);
00466 
00467   /* ----- Conclusion ------------------------------ */
00468 
00469   rms = sqrt(rms/npts);
00470   if (DEBUG) print_msg("Fit of correlation peak: status=%d, RMS=%g",
00471                        status,rms);
00472 
00473   if ((par[0] < parlb[0]+TOO_CLOSE) || (par[0] > parub[0]-TOO_CLOSE) || 
00474       (par[1] < parlb[1]+TOO_CLOSE) || (par[1] > parub[1]-TOO_CLOSE)) {
00475     /* We are within TOO_CLOSE px of the edge: something went wrong */
00476     print_warning("Fitting of X-peak dubious (peak closer than %.1f px from border)",
00477                   TOO_CLOSE);
00478   }
00479 
00480   /* Update parameters */
00481   par[0] -= winx;
00482   par[1] -= winy;                            /* Position offset wrt center */
00483   par[2] *= Imax; epar[2] *= Imax;           /* Intensity */
00484   par[5] *= Imax; epar[5] *= Imax;           /* Bkgnd */
00485 
00486   if (DEBUG && !status) {
00487     npix[0] = 2*winx;
00488     npix[1] = 2*winy;
00489 
00490     /* Dump central peak (use corr as storage place) */
00491     for (k=0, i=0; i<npix[0]; i++) 
00492       for (j=0; j<npix[1]; j++, k++) 
00493         corr[i][j] = y[k];
00494 
00495     sprintf(dbgname, "%s_peak%0d", DBG_PREFIX,index);
00496     print_msg("DEBUG: Storing arc X-peak (%dx%d) in %s",
00497               npix[0],npix[1],dbgname);
00498     disable_erase_flag();
00499     dump_2Darray((void **)corr,DOUBLE,npix,dbgname);
00500     restore_erase_flag();
00501 
00502     /* Dump central peak fit (use corr as storage) */
00503     for (k=0, i=0; i<npix[0]; i++) 
00504       for (j=0; j<npix[1]; j++, k++) 
00505         corr[i][j] = yf[k];
00506 
00507     sprintf(dbgname, "%s_peakfit%0d", DBG_PREFIX,index);
00508     print_msg("DEBUG: Storing arc X-peak fit (%dx%d) in %s",
00509               npix[0],npix[1],dbgname);
00510     disable_erase_flag();
00511     dump_2Darray((void **)corr,DOUBLE,npix,dbgname);
00512     restore_erase_flag();
00513   }
00514 
00515   /* Memory management */
00516   free2d(&corr,DOUBLE);
00517   free(x); free(y); free(yf);
00518 
00519   index++;
00520   
00521   return(status);
00522 }
00523 

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