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

extract/lib/proc_mpup.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 #ifdef HAVE_LIBDISLIN
00021 #include <dislin.h>
00022 #endif
00023 
00025 #define DISLIN 1
00026 
00028 #define MERGENEIGHBOR 1
00029 
00031 #define MEDSIG2GAUSSIG 1.4826
00032 
00033 /* === Doxygen Comment ======================================= */
00044 /* =========================================================== */
00045 
00046 int discard_null_mpups(int npup, int **ilist, int **jlist, float **value,
00047                        int ppup[])
00048 {
00049   int i,j, new_npup=npup;
00050 
00051   for (i=0; i<new_npup; i++) {
00052     if (ppup[i] == 0) { /* Discard this class: bring higher classes down 1 */
00053       new_npup --;
00054       for (j=i; j<new_npup; j++) {
00055         memcpy((int *)ilist[j], (int *)ilist[j+1], ppup[j+1]*sizeof(int));
00056         memcpy((int *)jlist[j], (int *)jlist[j+1], ppup[j+1]*sizeof(int));
00057         memcpy((double *)value[j], (double *)value[j+1], 
00058                ppup[j+1]*sizeof(double));
00059         ppup[j]  = ppup[j+1];
00060       }
00061       i--;              /* Warning: next class should be tested too! */
00062     }
00063   }
00064   /* We could reallocate arrays, but there's no realloc2d... */
00065 
00066   return(new_npup);
00067 }
00068 
00069 /* === Doxygen Comment ======================================= */
00081 /* =========================================================== */
00082 
00083 int detect_mpups(IMAGE2D mpup, double threshold, 
00084                  int *ppup[], double *fpup[], float *xpup[], float *ypup[])
00085 {
00086   int npup,i,j,k,l;
00087   int maxnpup, npixpup, flag;
00088   int minppup, maxppup;
00089   float mindist, maxdist, intens, di, dj, ftmp;
00090   int **ilist, **jlist, *work2;
00091   float **value;
00092   double *dtmp, *work1;
00093   double med,sigmed;
00094 
00095   maxnpup = 100*SNIFS_LENS_NX*SNIFS_LENS_NY;
00096   npixpup = 5000;                /* More than a full dead column */
00097   mindist  = 0.66*SNIFS_LENS_PX; /* 2/3 of the interlens distance */
00098   maxdist  = 2*mindist;          /* 4/3 of the interlens distance */
00099 
00100   alloc2d(&ilist,maxnpup,npixpup,INT);       /* coord. of pixels in mpup */
00101   alloc2d(&jlist,maxnpup,npixpup,INT);
00102   alloc2d(&value,maxnpup,npixpup,FLOAT);     /* value of pixels in mpup  */
00103 
00104   *ppup = (int *)malloc(maxnpup*sizeof(int)); /* nb of pixels in mpup     */
00105 
00106   /* ===== Detection ============================== */
00107 
00108   reset_print_progress();
00109   for (npup=0, j=0; j<mpup.ny; j++) { /* loop on rows first (natural order) */
00110     if (!DEBUG) print_progress("Finding classes", 100*(j+1)/mpup.ny, 5.);
00111     for (i=0; i<mpup.nx; i++) {                          /* loop on columns */
00112       /* Loop through all the mpup image */
00113       if ((intens = RD_frame(&mpup,i,j)) >= threshold) { /* something there */
00114         /* There's something here: does it belong to a previous mpup? */
00115         for (flag=0, k=npup-1; k>=0 && !flag; k--) {
00116           if ((ftmp = hypot((ilist[k][0]-i),(jlist[k][0]-j))) <= mindist) {
00117             /* It could belong to a previous mpup */
00118             for (l=0; l<(*ppup)[k] && !flag; l++) {
00119               if ((ABS(ilist[k][l] - i) <= MERGENEIGHBOR) && 
00120                   (ABS(jlist[k][l] - j) <= MERGENEIGHBOR)) {
00121                 /* It does belong to a previous mpup: update it */
00122                 flag = 1;
00123                 ilist[k][(*ppup)[k]] = i;
00124                 jlist[k][(*ppup)[k]] = j;
00125                 value[k][(*ppup)[k]] = intens;
00126                 (*ppup)[k]++;
00127 
00128                 if ((*ppup)[k] >= npixpup) {
00129                   print_error("There's more mpup pixels than expected (%d)",
00130                               npixpup);
00131                   exit_session(UNKNOWN);
00132                 }
00133               }
00134             }
00135           }
00136         }
00137         if (!flag) {                         /* new mpup */
00138           /* It doesn't belong to a previous mpup: create a new one */
00139           ilist[npup][0] = i;
00140           jlist[npup][0] = j;
00141           value[npup][0] = intens;
00142           (*ppup)[npup]  = 1;
00143           npup++;                            /* nb of mpup up to now */
00144 
00145           if (npup >= maxnpup) {
00146             print_error("There's more mpups than expected (%d)",maxnpup);
00147             exit_session(UNKNOWN);
00148           }
00149         } /* new mpup        */
00150       }   /* something there */
00151     }     /* loop on columns */
00152   }       /* loop on rows    */
00153   print_msg(" -> %d classes detected",npup);
00154 
00155   /* ===== Merging ============================== */
00156 
00157   print_msg("Merging classes closer than %dx%d px...",
00158             3*MERGENEIGHBOR,3*MERGENEIGHBOR);
00159   for (j=0; j<npup; j++) {                   /* Loop on all mpups */
00160     for (k=j+1; k<npup; k++) {               /* Loop on non-merged mpups */
00161       if ((ftmp = hypot((ilist[k][0]-ilist[j][0]),(jlist[k][0]-jlist[j][0]))) 
00162           <= mindist) {
00163         /* It could belong to a previous mpup */
00164         for (flag=0, i=0; i<(*ppup)[j] && !flag; i++) {
00165           for (l=0; l<(*ppup)[k] && !flag ; l++) {
00166             if ((ABS(ilist[k][l] - ilist[j][i]) <= 3*MERGENEIGHBOR) &&
00167                 (ABS(jlist[k][l] - jlist[j][i]) <= 3*MERGENEIGHBOR)) 
00168               flag = 1;                      /* Let's merge j and k! */
00169           } 
00170         } 
00171         if (flag) {                          /* Merge classe k into j */
00172 #ifdef MEGAVERBOSE
00173           if (DEBUG)
00174             print_msg("Classe %d (%dx%d, %d px) merged to "
00175                       "classe %d (%dx%d, %d px)",
00176                       k+1,ilist[k][0],jlist[k][0],(*ppup)[k],
00177                       j+1,ilist[j][0],jlist[j][0],(*ppup)[j]);
00178 #endif
00179           for (l=0; l<(*ppup)[k]; l++) {     /* Update pixels of j */
00180             ilist[j][(*ppup)[j]+l] = ilist[k][l];
00181             jlist[j][(*ppup)[j]+l] = jlist[k][l];
00182             value[j][(*ppup)[j]+l] = value[k][l];
00183           }
00184           (*ppup)[j] += (*ppup)[k];          /* Update nb of pixels in j*/
00185           (*ppup)[k] = 0;                    /* Discard merged classe k */
00186         }
00187       }
00188     } /* k-loop on non-merged classes */
00189   }   /* j-loop on all classes */
00190 
00191   npup = discard_null_mpups(npup, ilist, jlist, value, *ppup);
00192   print_msg(" -> %d classes after merging",npup);
00193 
00194   /* ===== Discarding strange classes ============================== */
00195 
00196   dtmp  = (double *)malloc(npup*sizeof(double));
00197   work1 = (double *)malloc(npup*sizeof(double));
00198   work2 =    (int *)malloc(npup*sizeof(int));
00199   for (i=0; i<npup; i++) dtmp[i] = (*ppup)[i];
00200 
00201   med = median_sig(dtmp, npup, &sigmed, work1, work2);
00202 
00203   free(dtmp); free(work1); free(work2);
00204 
00205   print_msg("Classe size: median = %.1f px, sigma = %.1f px",
00206             med,sigmed);
00207   sigmed *= 5*MEDSIG2GAUSSIG;
00208   minppup = MAX(0,(int)med - sigmed);
00209   maxppup =       (int)med + sigmed;
00210   print_msg("Discarding classes of size out of [%d,%d] px",minppup,maxppup);
00211 
00212   for (k=0; k<npup; k++) {
00213     if ((*ppup)[k] < minppup || (*ppup)[k] > maxppup) {
00214       if (DEBUG) 
00215         print_msg("Classe %d (%dx%d, %d px) size-discarded",
00216                   k+1,ilist[k][0],jlist[k][0],(*ppup)[k]);
00217       (*ppup)[k] = 0;
00218     }
00219   }
00220 
00221   npup = discard_null_mpups(npup, ilist, jlist, value, *ppup);
00222   print_msg(" -> %d classes after size selection",npup);
00223 
00224   *ppup =   (int *)realloc(*ppup,npup*sizeof(int)); /* nb of pixels */
00225   *fpup =(double *)malloc(npup*sizeof(double));     /* flux */
00226   *xpup = (float *)malloc(npup*sizeof(float));      /* centroid */
00227   *ypup = (float *)malloc(npup*sizeof(float));
00228 
00229   reset_print_progress();
00230   for (i=0; i<npup; i++) {     /* loop over mpups */
00231     if (!DEBUG) print_progress("Computing barycenters",100*(i+1)/npup, 5.);
00232     for (intens=di=dj=0, j=0; j<(*ppup)[i]; j++) { /* loop over pix in mpup */
00233       intens += value[i][j];
00234       di     += ilist[i][j]*value[i][j];
00235       dj     += jlist[i][j]*value[i][j];
00236     }
00237     (*fpup)[i] = intens;
00238     (*xpup)[i] = mpup.startx + (di/intens)*mpup.stepx;
00239     (*ypup)[i] = mpup.starty + (dj/intens)*mpup.stepy;
00240   }
00241   free2d(&ilist,INT); free2d(&jlist,INT); free2d(&value,FLOAT);
00242 
00243   return(npup);
00244 }
00245 
00246 /* === Doxygen Comment ======================================= */
00276 /* =========================================================== */
00277 
00278 int study_triangulation(IMAGE2D mpup, float xpup[], float ypup[], int npup,
00279                         int iadj[], int iend[], double xtri[], double ytri[], 
00280                         double stri[], double ltri[], double atri[])
00281 {
00282   int *lnb;
00283   
00284   int i,j,k,i1,i2,ntri,nside,npts, flag;
00285   float mindist, maxdist, maxhp;
00286   float x0,y0,x1,y1,x2,y2,dx,dy;
00287   float dist01, dist02, dist12, halfper;
00288   
00289   mindist = 0.66*SNIFS_LENS_PX; /* 2/3 of the interlens distance */
00290   maxdist = 2*mindist;          /* 4/3 of the interlens distance */
00291   maxhp   = 1.85*SNIFS_LENS_PX; /* between 1+sqrt(2)/2 and 2 (see below) */
00292 
00293 #ifdef HAVE_LIBDISLIN 
00294   if (DEBUG && !DISLIN) 
00295     print_warning("Triangulation graphical output desactivated (%s:%d)",
00296                   __FILE__,__LINE__);
00297   if (DEBUG && DISLIN) {
00298     Plot plot;
00299     char title[4*lg_ident+1];
00300     /* Graphic output initialization */
00301     sprintf(title,"dislin_triangulation.eps");
00302     print_msg("DISLIN: Output graphics in %s",title);
00303     setfil(title);                           /* Output filename */
00304 
00305     plot_initialize(&plot, 1, "EPS", NULL, &(plot.Symbol), 0,0);
00306     plot_set_limits(xpup,ypup,npup,&(plot.Limits));
00307     plot_increase_limits(&(plot.Limits));
00308     plot_start(&plot,0);
00309     plot_axes(&plot,"x [px]","y [px]");
00310 
00311     sprintf(title,"1st-order mpup triangulation\n\nFrame %s, %d nodes",
00312             mpup.name,npup);
00313     plot_title(title);
00314   }
00315 #endif /* HAVE_LIBDISLIN */
00316 
00317   reset_print_progress();
00318   for (ntri=0,nside=0,lnb=iadj,i=0; i<npup; i++) { /* loop on nodes */
00319     print_progress("Compute triangle surface & orientation",100*(i+1)/npup,5.);
00320 
00321     /* Number of neighbours around node i */
00322     if (i) npts = iend[i]-iend[i-1];
00323     else   npts = iend[0];
00324 
00325     for (j=0; j<npts; j++) { /* Loop over the triangles around node i */
00326       i1 = lnb[j];
00327       if (j != npts-1) i2 = lnb[j+1]; /* take next neighbor as 3rd point */
00328       else             i2 = lnb[0];   /* loop is over: take 1st neighbor */
00329       if ((i1 == 0) || (i2 == 0)) {   /* node i is at the boundary */
00330 #ifdef HAVE_LIBDISLIN
00331         /* plot outer mpup */
00332         if (DEBUG && DISLIN) { 
00333           color("red"); 
00334           rlsymb(1,xpup[i],ypup[i]); 
00335         }
00336 #endif
00337         continue; 
00338       }
00339 
00340       /* Corners of the triangle */
00341       x0 = xpup[i]; x1 = xpup[i1-1]; x2 = xpup[i2-1];
00342       y0 = ypup[i]; y1 = ypup[i1-1]; y2 = ypup[i2-1];
00343 
00344       /* Side lengthes and 1/2-perimeter */
00345       dist01 = hypot((x0-x1),(y0-y1));
00346       dist02 = hypot((x0-x2),(y0-y2));
00347       dist12 = hypot((x2-x1),(y2-y1));
00348       halfper = (dist01 + dist02 + dist12)/2.;
00349 
00350       /* Discriminating the "good" triangles (that is corresponding to
00351        half/squares): 1- *exactly* 2 sides of the triangle should be
00352        of order of the interlens distance, 2- halfper should be of
00353        order 1+sqrt(2)/2 ~ 1.7 the interlens distance to remove ~flat
00354        triangles. */
00355       flag = 0;
00356       if ((dist01 < mindist) || (dist01 > maxdist)) flag ++;
00357       if ((dist02 < mindist) || (dist02 > maxdist)) flag ++;
00358       if ((dist12 < mindist) || (dist12 > maxdist)) flag ++;
00359       if ((flag != 1) || (halfper > maxhp)) continue;
00360 
00361       /* Triangle barycenter */
00362       xtri[ntri] = (x0 + x1 + x2)/3;
00363       ytri[ntri] = (y0 + y1 + y2)/3;
00364 
00365       /* Discard previously found triangles */
00366       for (flag=0,k=0; k<ntri; k++) {
00367         if (hypot(xtri[ntri]-xtri[k],ytri[ntri]-ytri[k]) < 1e-10) {
00368           flag = 1;
00369           break;
00370         }
00371       }
00372       if (flag) continue;
00373 
00374 #ifdef HAVE_LIBDISLIN
00375       /* draw triangle */
00376       if (DEBUG && DISLIN) {
00377         color("blue"); 
00378         rline(x0,y0,x1,y1);rline(x0,y0,x2,y2);rline(x1,y1,x2,y2);
00379       }
00380 #endif
00381 
00382       /* Triangle surface: with side length a,b,c, and p = (a+b+c)/2,
00383          S = sqrt( p*(p-a)*(p-b)*(p-c) ) */
00384       stri[ntri] = sqrt(halfper * 
00385                         (halfper-dist01)*(halfper-dist02)*(halfper-dist12));
00386 
00387       /* Triangle side length and orientation */
00388       if ((dist01 > dist02) && (dist01 > dist12)) { /*   0 */
00389         dx = x2 - x1;                               /*  /| */
00390         dy = y2 - y1;                               /* 1-2 */
00391         ltri[nside++] = dist02;
00392         ltri[nside++] = dist12;
00393       }
00394       if ((dist02 > dist01) && (dist02 > dist12)) { /*   2 */
00395         dx = x1 - x0;                               /*  /| */
00396         dy = y1 - y0;                               /* 0-1 */
00397         ltri[nside++] = dist01;
00398         ltri[nside++] = dist12;
00399       }
00400       if ((dist12 > dist01) && (dist12 > dist02)) { /*   2 */
00401         dx = x0 - x1;                               /*  /| */
00402         dy = y0 - y1;                               /* 1-0 */
00403         ltri[nside++] = dist01;
00404         ltri[nside++] = dist02;
00405       }
00406       if (ABS(dx) >= ABS(dy)) atri[ntri] = atan(dy/dx);
00407       else                    atri[ntri] = atan(-dx/dy);
00408 
00409       ntri++; /* one more valid triangle */
00410     } /* Loop over the triangles around node i */
00411     lnb += npts; /* next set of local neighbors */
00412 
00413 #ifdef HAVE_LIBDISLIN
00414     /* Number mpup */
00415     if (DEBUG && DISLIN) { 
00416       color("fore"); 
00417       rlnumb(i+1,-1,xpup[i],ypup[i]); 
00418     }
00419 #endif
00420   } /* loop on nodes */
00421 
00422 #ifdef HAVE_LIBDISLIN
00423   /* Graphic output conclusion */
00424   if (DEBUG && DISLIN) disfin();
00425 #endif
00426 
00427   return(ntri);  
00428 }
00429 
00430 /* === Doxygen Comment ======================================= */
00440 /* =========================================================== */
00441 
00442 void mpup_orders(double fpup[], int opup[], int npup, 
00443                  double f[3], double df[3])
00444 {
00445   int i, n0,n1,n2;
00446   double *f0,*f1,*f2;
00447   
00448   f0 = (double *)malloc(npup*sizeof(double));
00449   f1 = (double *)malloc(npup*sizeof(double));
00450   f2 = (double *)malloc(npup*sizeof(double));
00451   for (n0=n1=n2=0,i=0; i<npup; i++) {
00452     switch(opup[i]) {
00453     case 0: f0[n0++] = fpup[i]; break;
00454     case 1: f1[n1++] = fpup[i]; break;
00455     case 2: f2[n2++] = fpup[i]; break;
00456     }
00457   }
00458   mean_sig(f0,n0,f,  df);
00459   mean_sig(f1,n1,f+1,df+1);
00460   mean_sig(f2,n2,f+2,df+2);
00461   print_msg("   %d 1st order mpups (flux = %.2f +- %.2f)\n"
00462             "   %d 0th-order mpups (flux = %.2f +- %.2f)\n"
00463             "   %d 2nd-order mpups (flux = %.2f +- %.2f)",
00464             n1,f[1],df[1],n0,f[0],df[0],n2,f[2],df[2]);
00465   free(f0); free(f1); free(f2);
00466 }
00467 
00468 /* === Doxygen Comment ======================================= */
00481 /* =========================================================== */
00482 
00483 void rotate_coord(float x[],float y[],int n,
00484                   float x0,float y0,float angle, float xr[],float yr[])
00485 {
00486   int i;
00487   float t;
00488   double c=cos(angle),s=sin(angle);
00489   
00490   for (i=0; i<n; i++) {
00491     xr[i] = x[i] - x0; /* centering */
00492     yr[i] = y[i] - y0;
00493     t = xr[i];         /* rotating */
00494     xr[i] = xr[i]*c - yr[i]*s;
00495     yr[i] =     t*s + yr[i]*c;
00496     xr[i] += x0;       /* offseting */
00497     yr[i] += y0;
00498   }
00499 }
00500 
00501 /* === Doxygen Comment ======================================= */
00515 /* =========================================================== */
00516 
00517 int select_mpup_dist(float x0, float y0, float maxdist, 
00518                      int npup, float xpup[], float ypup[], int opup[],
00519                      int icen[]) 
00520 {
00521   int ncen, i;
00522   float ftmp;
00523 
00524   print_msg("Mpup selection: 1st-order mpups %s to (%.2f,%.2f) "
00525             "than %.2f px",
00526             (maxdist > 0) ? "closer" : "further",x0,y0,ABS(maxdist));
00527 
00528   for (ncen=i=0; i<npup; i++) {
00529     if (opup[i] != 1) continue;      /* Keep only 1st-order mpups */
00530     ftmp = hypot(xpup[i] - x0,ypup[i] - y0);
00531     if ((maxdist>0 && ftmp<ABS(maxdist)) ||
00532         (maxdist<0 && ftmp>ABS(maxdist))) {
00533       icen[ncen++] = i;
00534     }
00535   }
00536   return(ncen);
00537 }
00538 
00539 /* === Doxygen Comment ======================================= */
00552 /* =========================================================== */
00553 
00554 int select_mpup_sharpest_smooth(float maxdist,
00555                                 int npup,float xpup[],float ypup[], 
00556                                 int opup[], float spup[],
00557                                 int icen[]) 
00558 {
00559   int *ipup1, *hpup1;
00560   int ncen, i,k, npup1, imin,imax;
00561 
00562   float *xpup1,*ypup1,*spup1;
00563   float smin,smax;
00564 
00565   /* Selection of the 1st-order mpups */
00566   xpup1 = (float *)malloc(npup*sizeof(float));
00567   ypup1 = (float *)malloc(npup*sizeof(float)); /* Coordinates */
00568   spup1 = (float *)malloc(npup*sizeof(float)); /* Dispersion */
00569   ipup1 =   (int *)malloc(npup*sizeof(int));   /* Index */
00570   hpup1 =   (int *)malloc(npup*sizeof(int));   /* Nb of neighbors */
00571   for (npup1=i=0; i<npup; i++) {
00572     if (opup[i] != 1) continue;
00573     xpup1[npup1] = xpup[i];
00574     ypup1[npup1] = ypup[i];
00575     spup1[npup1] = spup[i];
00576     ipup1[npup1] = i;
00577     npup1++;
00578   }
00579   print_msg("Mpup selection: filtering %d 1st-order mpups",npup1);
00580 
00581   /* Spatial filtering */
00582   spatial_filter(npup1,xpup1,ypup1,spup1, maxdist,'M',4,hpup1);
00583 
00584   /* Keeping mpups w/ full set of neighbors */
00585   for (k=i=0; i<npup1; i++) {
00586     if (hpup1[i] < 9) continue;
00587     xpup1[k] = xpup1[i];
00588     ypup1[k] = ypup1[i];
00589     spup1[k] = spup1[i];
00590     ipup1[k] = ipup1[i];
00591     k++;
00592   }
00593   npup1 = k;
00594   print_msg("Mpup selection: %d filtered 1st-order mpups",npup1);
00595 
00596   /* Min/max filtered dispersions */
00597   minmax_f(spup1,npup1, &smin,&smax, &imin,&imax);
00598   imin = ipup1[imin];
00599   imax = ipup1[imax];
00600   free(xpup1); free(ypup1); free(spup1); free(ipup1); free(hpup1);
00601 
00602   print_msg("Mpup selection:\n"
00603             "   min. filtered X-disp.: %.2f at (%.2fx%.2f)\n"
00604             "   max. filtered X-disp.: %.2f at (%.2fx%.2f)",
00605             smin,xpup[imin],ypup[imin],smax,xpup[imax],ypup[imax]);
00606 
00607   /* Looking for neighbors */
00608   ncen = select_mpup_dist(xpup[imin],ypup[imin],maxdist, 
00609                           npup,xpup,ypup,opup, icen);
00610 
00611   return(ncen);
00612 }

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