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

lib/stats.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 <gsl/gsl_histogram.h>
00021 
00022 /* === Doxygen Comment ======================================= */
00038 /* =========================================================== */
00039 
00040 void frame_stats(const IMAGE2D *frame, double *mean, double *sigma, 
00041                  double *median, double *sigmed) 
00042 {
00043   int *iwork;
00044   int i,j,npts;
00045   double dtmp, sum, sum2;
00046   double *vec, *work;
00047 
00048   /* Allocation */
00049   npts = (frame->nx)*(frame->ny);
00050   vec  = (double *)malloc(npts*sizeof(double));
00051 
00052   /* Preparation */
00053   sum  = 0.;
00054   sum2 = 0.;
00055   for (j=0; j<frame->ny; j++) { /* Image-oriented loop (i running faster) */
00056     for (i=0; i<frame->nx; i++) {
00057       dtmp  = (double)RD_frame(frame,i,j);
00058       sum  += dtmp;
00059       sum2 += dtmp*dtmp;
00060       vec[i*(frame->ny) + j] = dtmp;
00061     }
00062   }
00063 
00064   /* Mean & sigma (cf. p.58, Elements of Math. Stat.) */
00065   *mean  = sum/npts;
00066   *sigma = sqrt(sum2/(npts-1) - (*mean)*(sum/(npts-1)));
00067 
00068   /* Median & median sigma */
00069   if (median!=NULL && sigmed!=NULL) {
00070     work = (double *)malloc(npts*sizeof(double));
00071     iwork=    (int *)malloc(npts*sizeof(int));
00072 
00073     *median = median_sig(vec,npts,sigmed,work,iwork);
00074 
00075     free(work); free(iwork);
00076   }
00077 
00078   free(vec); 
00079 }
00080 
00081 /* === Doxygen Comment ======================================= */
00097 /* =========================================================== */
00098 
00099 SPECTRUM *array_histo(const double *array, int npts, 
00100                       double min, double max, int nbin)  
00101 {
00102   gsl_histogram *histo;
00103   SPECTRUM *histospec;
00104   int nb,i;
00105   double start,step;
00106   double *range;
00107 
00108   /* Ranges of the histogram (one more than bins) */
00109   step = (max-min)/ABS(nbin);
00110   if (nbin<=0) {
00111     /* If nbin<0, total nbin = ABS(nbin) + 2 */
00112     nb = ABS(nbin)+2;
00113     range = (double *)malloc((nb+1)*sizeof(double));
00114     range[0] = -MAXDOUBLE;
00115     for (i=1; i<nb; i++) range[i] = min+(i-1)*step;
00116     range[nb] = MAXDOUBLE;
00117     start = min - step/2.;
00118   }
00119   else {
00120     /* If nbin>0, total nbin = nbin */
00121     nb = nbin;
00122     range = (double *)malloc((nb+1)*sizeof(double));
00123     for (i=0; i<=nb; i++) range[i] = min+i*step;
00124     start = min + step/2.;
00125   }
00126   
00127   /* Histogram initialization */
00128   histo = gsl_histogram_alloc(nb);
00129   gsl_histogram_set_ranges(histo, range, nb+1);
00130 
00131   /* Histogram construction */
00132   for (i=0; i<npts; i++) gsl_histogram_increment(histo,array[i]);
00133 
00134   /* Store the histogram in a spectrum */
00135   histospec = (SPECTRUM *)malloc(sizeof(SPECTRUM));
00136   create_spec_mem(histospec, nb, start, step, FLOAT);
00137   for (i=0; i<nb; i++) WR_spec(histospec,i,(float)gsl_histogram_get(histo,i));
00138 
00139   /* Memory management */
00140   gsl_histogram_free(histo);
00141   free(range);
00142 
00143   return(histospec);
00144 }
00145 
00146 /* === Doxygen Comment ======================================= */
00166 /* =========================================================== */
00167 
00168 SPECTRUM *subframe_histo(const IMAGE2D *frame, 
00169                          int imin, int imax, int jmin, int jmax, 
00170                          double min, double max, int nbin)  
00171 {
00172   SPECTRUM *histospec;
00173   int npts,i,j,k;
00174   double *array;
00175 
00176   npts = (imax-imin+1)*(jmax-jmin+1);
00177 
00178   /* Copy frame into array */
00179   array = (double *)malloc(npts*sizeof(double));
00180   for (k=0, j=jmin; j<=jmax; j++) /* Image-oriented loop (i running faster) */
00181     for (i=imin; i<=imax; i++, k++) 
00182       array[k] = (double)RD_frame(frame,i,j);
00183 
00184   /* Array histogram */
00185   histospec = array_histo(array,npts,min,max,nbin);
00186 
00187   free(array);
00188 
00189   return(histospec);
00190 }
00191 
00192 /* === Doxygen Comment ======================================= */
00205 /* =========================================================== */
00206 
00207 SPECTRUM *frame_histo(const IMAGE2D *frame, double min, double max, int nbin)  
00208 {
00209   SPECTRUM *histospec;
00210 
00211   histospec = subframe_histo(frame,0,frame->nx-1,0,frame->ny-1, min,max,nbin);
00212 
00213   return(histospec);
00214 }
00215 
00216 /* === Doxygen Comment ======================================= */
00227 /* =========================================================== */
00228 
00229 double histo_threshold(const SPECTRUM *histospec, float frac)
00230 {
00231   int i;
00232   double sum,limit,dx,threshold;
00233   
00234   /* Total sum of the histogram */
00235   for (i=0, sum=0; i<histospec->npts; i++) sum += (float)RD_spec(histospec,i);
00236 
00237   /* Compute the threshold */
00238   limit = frac*sum;
00239   for (i=0, sum=0; i<histospec->npts; i++) {
00240     sum += RD_spec(histospec,i);
00241     if (sum>=limit) {
00242       /* Linear interpolation over the last histogram bin */
00243       dx = (sum-limit)/RD_spec(histospec,i);
00244       threshold = histospec->start + (i+0.5-dx)*histospec->step;
00245       break;
00246     }
00247   }
00248   
00249   return(threshold);
00250 }
00251 
00252 /* === Doxygen Comment ======================================= */
00260 /* =========================================================== */
00261 
00262 float median_f(const float x[], int n)
00263 {
00264   int *warr,i;
00265   float med;
00266   double *xd;
00267   
00268   warr =  (int *)malloc(n*sizeof(int));
00269   xd = (double *)malloc(n*sizeof(double));
00270 
00271   for (i=0; i<n; i++) xd[i] = x[i];
00272 
00273   med = median(xd,n,warr);
00274   free(warr); free(xd); 
00275 
00276   return(med);
00277 }
00278 
00279 /* === Doxygen Comment ======================================= */
00291 /* =========================================================== */
00292 
00293 double mean_sig_f(const float x[], int n, double *mean, double *sigma)
00294 {
00295   int i;
00296   double sum,sum2;
00297   
00298   for (sum=0, sum2=0, i=0; i<n; i++) {
00299     sum  += x[i];
00300     sum2 += SQ(x[i]);
00301   }
00302   /* Mean & sigma (cf. p.58, Elements of Math. Stat.) */
00303   *mean  = sum/n;
00304   *sigma = sqrt(sum2/(n-1) - (*mean)*(sum/(n-1)));
00305 
00306   return(*mean);
00307 }
00308 
00309 /* === Doxygen Comment ======================================= */
00322 /* =========================================================== */
00323 
00324 void minmax(const double x[], int n, 
00325             double *min, double *max, int *imin, int *imax)
00326 {
00327   int i;
00328   
00329   for (*min=MAXDOUBLE, *max=-MAXDOUBLE, i=0; i<n; i++) {
00330     if (x[i] < *min) { *min = x[i]; *imin = i; }
00331     if (x[i] > *max) { *max = x[i]; *imax = i; }
00332   }
00333 }
00334 
00335 /* === Doxygen Comment ======================================= */
00348 /* =========================================================== */
00349 
00350 void minmax_f(const float x[], int n, 
00351               float *min, float *max, int *imin, int *imax)
00352 {
00353   int i;
00354   
00355   for (*min=MAXFLOAT, *max=-MAXFLOAT, i=0; i<n; i++) {
00356     if (x[i] < *min) { *min = x[i]; *imin = i; }
00357     if (x[i] > *max) { *max = x[i]; *imax = i; }
00358   }
00359 }
00360 
00361 /* === Doxygen Comment ======================================= */
00372 /* =========================================================== */
00373 
00374 int locate_sorted(const double x[], int n, double x0)
00375 {
00376   int il,iu,im,order;
00377 
00378   il = -1;
00379   iu = n;
00380 
00381   order = (x[n-1] >= x[0]); /* Ascending or descending order? */
00382 
00383   while (iu-il > 1) {
00384     im = (il+iu)/2;
00385     if ((x0 >= x[im]) == order) il = im;
00386     else iu = im;
00387   }
00388   
00389   return(il);
00390 }
00391 
00392 /* === Doxygen Comment ======================================= */
00404 /* =========================================================== */
00405 
00406 void invert_indexx(int idx[], int n)
00407 {
00408   int *xdi;
00409   int i;
00410   
00411   xdi = (int *)malloc(n*sizeof(int));
00412 
00413   for (i=0; i<n; i++) xdi[idx[i]] = i;
00414 
00415   memcpy((int *)idx,(int *)xdi,n*sizeof(int));
00416   free(xdi);
00417 }
00418 
00419 /* Replacement for Fit_polynom in IFU_C_mathlibs-6.1/C_libs/math/least_sq/misc_routines.c */
00420 
00421 /*-----------------------------------------------------------------------------
00422 !
00423 !.blk                          Polynomial curve fitting
00424 !
00425 !.func                            Fit_polynom()
00426 !
00427 !.purp              fit an orthogonal polynomials upon given intensities
00428 !.desc
00429 ! int Fit_polynom(x, y, nbp, input_degree, output_degree, coef, dy, khi2)
00430 !
00431 ! float  *x;           x coordinates
00432 ! float  *y;           y coordinates
00433 ! int    npts;         number of coordinates in previous arrays
00434 ! int    input_deg;    estimated degree of polynom
00435 ! int    *output_deg;  degree found for polynom (<= input_deg)
00436 ! float  *coef;        coefficients (array size = 3*(input_deg+1))
00437 ! float  *dy;          lsq residuals (allocated array size=npts)
00438 ! float  *khi2;        value of khi2
00439 !.ed
00440 ------------------------------------------------------------------------------*/
00441 
00442 /* === Doxygen Comment ======================================= */
00449 /* =========================================================== */
00450 
00451 int Fit_polynom2(float *x, float *y, int npts, int input_deg, int *output_deg,
00452                  float *coef, float *dy, float *chi2)
00453 {
00454   int const max_deg=30;
00455 
00456   int iff, i, i1, n, k, imax, END;
00457 
00458   double w[max_deg], p[max_deg];
00459   double *a,*b,*s, *dcoef;
00460   double xnu, delta, p2, dchi, f, f95;
00461 
00462   /*   ------------------------------------------------------------------------
00463        Least square fit of data set (X,Y) to orthogonal polynomials,
00464        Y(X) = SUM[ S(I)*P(I,X) ],
00465        of degree OUTPUT_DEG<30, which are recovered by recursion:
00466        P(1,X)   = 1
00467        P(2,X)   = X-A(1)
00468        P(N+1,X) = (X-A(N))*P(N,X)-B(N)*P(N-1,X)
00469 
00470        On entrance OUTPUT_DEG sets the maximum order to the fit, a lower order
00471        is returned if two consecutive higher order coefficients S(I) fail
00472        the 95 % significance of the F test.
00473        On exit DY contains the residual vector, and CHI2 is the summed squared
00474        residuals
00475        -------------------------------------------------------------------------- */
00476 
00477   if (input_deg+1 > max_deg) return(1);
00478   imax = MAX(input_deg+1,2);
00479 
00480   memcpy((char *)dy,(char *)y,npts*sizeof(float)); /* set dy[] with y[] val */
00481   dcoef = (double *)malloc(3*(input_deg+1)*sizeof(double)); /* Work in double */
00482 
00483   a = dcoef;
00484   b = dcoef +   (input_deg+1);
00485   s = dcoef + 2*(input_deg+1);
00486   *output_deg = input_deg;
00487 
00488   /* Initialization */
00489   for (i=0; i<imax; i++) a[i] = b[i] = s[i] = w[i] = 0;
00490 
00491   /* zero order */
00492   i = 0;
00493   p[0] = 1.;
00494   for (n=0; n<npts; n++) {
00495     w[0] += 1.;
00496     s[0] += dy[n];
00497     a[0] += x[n];
00498   }
00499   s[0] /= w[0];
00500   a[0] /= w[0];
00501 
00502   /* Higher orders */
00503   xnu = npts-1;
00504   END = 0;
00505 
00506   while (!END) {
00507 
00508     iff = 1;
00509     while (iff <= 2) {
00510       i1 = i;
00511       if (i+1 < imax) i++;
00512       *chi2 = 0;
00513       for (n=0; n<npts; n++) {
00514         p[1] = x[n]-a[0];
00515         if (i!=1) {
00516           for (k=2; k<=i; k++) {
00517             p[k] = (x[n]-a[k-1])*p[k-1] - b[k-1]*p[k-2];
00518           }
00519         }
00520         delta = dy[n] - s[i1]*p[i1];
00521         dy[n] = delta;
00522         *chi2 += SQ(delta);
00523         if (imax <= i1+1) continue;
00524         p2 = SQ(p[i]);
00525         s[i] += delta*p[i];
00526         a[i] += x[n]*p2;
00527         w[i] += p2;
00528       }
00529       if (imax <= i1+1) {
00530         END = 1; break;
00531       }
00532       a[i] /= w[i];
00533       b[i] = w[i]/w[i1];
00534       s[i] /= w[i];
00535 
00536       /* F-test for significance */
00537       xnu -= 1.;
00538       dchi = SQ(s[i])*w[i];
00539       if (dchi >= *chi2) break;
00540       f = xnu*dchi/(*chi2-dchi);
00541       if (xnu == 0)
00542         f95 = 1.e30;
00543       else
00544         f95 = 3.84+(10.+(12.+(30.+105./xnu/xnu)/xnu)/xnu)/xnu;
00545       if (f > f95) break;
00546       xnu += 1.;
00547       s[i] = 0.;
00548       iff++;
00549     } /* while (iff <= 2) */
00550 
00551     if (iff > 2) END = 1;
00552 
00553   } /* while(!END) */
00554 
00555   *output_deg = MIN(imax-1,i1+1);
00556   *output_deg += -iff + 1;
00557 
00558   /* Copy back double array to float array */
00559   for (i=0; i<3*(input_deg+1); i++) coef[i] = dcoef[i];
00560 
00561   *chi2 = sqrt((double)(*chi2))/(npts-1);
00562   if (*output_deg < input_deg) {
00563     for (i=0; i<=*output_deg; i++) coef[(*output_deg+1)+i]   = b[i];
00564     for (i=0; i<=*output_deg; i++) coef[(*output_deg+1)*2+i] = s[i];
00565   }
00566 
00567   free(dcoef);
00568   
00569   return(0);
00570 }
00571 
00572 /* === Doxygen Comment ======================================= */
00581 /* =========================================================== */
00582 
00583 void filter_med_array(float array[], int npts, int radius, float f_array[])
00584 {
00585   int *work;
00586   int i,j,k,jmin,jmax,npw;
00587   double *window;
00588 
00589   npw = 2*radius+1;                          /* Filtering window full-width */
00590   window = (double*)malloc(npw*sizeof(double));
00591   work   =    (int*)malloc(npw*sizeof(int));
00592 
00593   for (i=0; i<npts; i++) {                   /* Median filtering */
00594     jmin = MAX(0,i-radius);
00595     jmax = MIN(npts-1, i+radius);            /* npts-1 because j<=jmax */
00596     for (k=0,j=jmin; j<=jmax; j++,k++) window[k] = array[j];
00597     f_array[i] = median(window,k,work);
00598   }
00599 
00600   free(window); free(work);
00601 }
00602 

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