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

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

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