imcore_background.c

00001 /*
00002 
00003 $Id: imcore_background.c,v 1.8 2007/12/19 13:16:50 jim Exp $
00004 
00005 */
00006 
00007 #include <stdio.h>
00008 #include <stdlib.h>
00009 #include <string.h>
00010 #include <math.h>
00011 
00012 #include <cpl.h>
00013 
00014 #include "floatmath.h"
00015 #include "util.h"
00016 #include "imcore.h"
00017 
00018 static int **hist = NULL;
00019 static int *nnp = NULL;
00020 static int npvx;
00021 static int npvy;
00022 static void tidy(void);
00023 static void sortit (float [], int);
00024 
00025 extern int imcore_background(ap_t *ap, int nbsize, float nullval) {
00026     float fracx,fracy,skymed,sigma,skymedc,sigmac,avsky,fnbsize,dely,delx;
00027     float t1,t2,dsky,*map,**bvals,*work;
00028     int ifracx,ifracy,nbsizx,nbsizy,nbx,nby,npixstripe,l,i,ll;
00029     int isquare,ilev,j,iclip,mcpix,iloop,irej,nbsizo2,kk,k,iby,ibyp1,ibx,ibxp1;
00030     int *shist,*conf;
00031     unsigned char *mflag;
00032     long nx,ny;
00033 
00034     /* Set up some variables */
00035 
00036     map = ap->indata;
00037     conf = ap->confdata;
00038     mflag = ap->mflag;
00039     nx = ap->lsiz;
00040     ny = ap->csiz;
00041 
00042     /* check to see if nbsize is close to exact divisor */
00043 
00044     fracx = ((float)nx)/((float)nbsize);
00045     fracy = ((float)ny)/((float)nbsize);
00046     ifracx = (int)(fracx + 0.1);
00047     ifracy = (int)(fracy + 0.1);
00048     nbsizx = nx/ifracx;
00049     nbsizy = ny/ifracy;
00050     nbsize = MAX(NINT(0.9*nbsize), MIN(nbsize, MIN(nbsizx,nbsizy)));
00051     nbsize = MIN(nx,MIN(ny,nbsize)); /* trap for small maps */
00052 
00053     /* Divide the map into partitions */
00054 
00055     nbx = nx/nbsize;
00056     nby = ny/nbsize;
00057     npixstripe = nbsize*nx;
00058     npvx = nbx;
00059     npvy = nby;
00060 
00061     /* Get histogram workspace if you can */
00062 
00063     hist = cpl_malloc(nbx*sizeof(int *));
00064     for (l = 0; l < nbx; l++) 
00065         hist[l] = cpl_malloc(MAXHIST*sizeof(int));
00066 
00067     /* Same for background values array */
00068 
00069     bvals = cpl_malloc(nby*sizeof(float *));
00070     for (l = 0; l < nby; l++) 
00071         bvals[l] = cpl_malloc(nbx*sizeof(float));
00072 
00073     /* Store some of this away for use later */
00074 
00075     ap->backmap.nbx = nbx;
00076     ap->backmap.nby = nby;
00077     ap->backmap.nbsize = nbsize;
00078     ap->backmap.bvals = bvals;
00079 
00080     /* Finally a counter array */
00081 
00082     nnp = cpl_malloc(nbx*sizeof(int));
00083 
00084     /* Loop for each row of background squares. Start by initialising
00085        the accumulators and histograms */
00086 
00087     for (l = 0; l < nby; l++) {
00088         memset((char *)nnp,0,nbx*sizeof(*nnp));
00089         for (i = 0; i < nbx; i++) 
00090             memset((char *)hist[i],0,MAXHIST*sizeof(int));
00091 
00092         /* Skim through the data in this stripe. Find out which square each
00093            belongs to and add it it to the relevant histogram */
00094 
00095         ll = l*npixstripe;
00096         for (i = 0; i < npixstripe; i++) {
00097             if (map[ll+i] != nullval && conf[ll+i] > 0 &&
00098                 mflag[ll+i] != MF_SATURATED) {
00099                 isquare = (int)((float)(i % nx)/(float)nbsize);
00100                 isquare = MIN(nbx-1,MAX(0,isquare));
00101                 ilev = MIN(MAXHISTVAL,MAX(MINHISTVAL,NINT(map[i+ll])));
00102                 hist[isquare][ilev-MINHISTVAL] += 1;
00103                 nnp[isquare] += 1;
00104             }
00105         }
00106   
00107         /* but only do background estimation if enough pixels ----------- */
00108 
00109         for (j = 0; j < nbx; j++) {
00110             if (nnp[j] > 0.25*nbsize*nbsize){
00111                 shist = hist[j];
00112                 imcore_medsig(shist,MAXHIST,MINHISTVAL-1,nnp[j],&skymed,&sigma);
00113 
00114                 /* do an iterative 3-sigma upper clip to give a more robust
00115                    estimator */
00116 
00117                 iclip = MAXHISTVAL;
00118                 mcpix = nnp[j];
00119                 skymedc = skymed;
00120                 sigmac = sigma;
00121                 for(iloop = 0; iloop < 3; iloop++) {
00122                     irej = 0;
00123                     for(i = NINT(skymedc+3.0*sigmac); i <= iclip; i++)
00124                         irej += shist[i-MINHISTVAL];
00125                     if (irej == 0)
00126                         break;
00127                     iclip = NINT(skymedc+3.0*sigmac) - 1;
00128                     mcpix = mcpix - irej;
00129                     imcore_medsig(shist,MAXHIST,MINHISTVAL-1,mcpix,&skymedc,
00130                                   &sigmac);
00131                 }
00132                 bvals[l][j] = skymedc;
00133             } else {
00134                 bvals[l][j] = -1000.0;
00135             }
00136         }
00137     }
00138 
00139     /* filter raw background values */
00140 
00141     bfilt(bvals,nbx,nby);
00142 
00143     /* compute average sky level */
00144 
00145     work = cpl_malloc(nbx*nby*sizeof(*work));
00146     k = 0;
00147     for(l = 0; l < nby; l++)
00148         for(j = 0; j < nbx; j++) 
00149             work[k++] = bvals[l][j];
00150     sortit(work,nbx*nby);
00151     avsky = work[(nbx*nby)/2];
00152     freespace(work);
00153 
00154     /* ok now correct map for background variations and put avsky back on */
00155 
00156     nbsizo2 = nbsize/2;
00157     fnbsize = 1.0/((float)nbsize);
00158     for (k = 0; k < ny; k++) {
00159         kk = k*nx;
00160 
00161         /* Nearest background pixel vertically */
00162 
00163         iby = (k + 1 + nbsizo2)/nbsize;
00164         ibyp1 = iby + 1;
00165         iby = MIN(nby,MAX(1,iby));
00166         ibyp1 = MIN(nby,ibyp1);
00167         dely = (k + 1 - nbsize*iby + nbsizo2)*fnbsize;
00168 
00169         for (j = 0; j < nx; j++) {
00170             if (map[kk+j] == nullval) 
00171                 continue;
00172 
00173             /* nearest background pixel across */
00174 
00175             ibx = (j + 1 + nbsizo2)/nbsize;
00176             ibxp1 = ibx + 1;
00177             ibx = MIN(nbx,MAX(1,ibx));
00178             ibxp1 = MIN(nbx,ibxp1);
00179             delx = (j + 1 - nbsize*ibx + nbsizo2)*fnbsize;
00180 
00181             /* bilinear interpolation to find background */
00182             
00183             t1 = (1.0 - dely)*bvals[iby-1][ibx-1] + dely*bvals[ibyp1-1][ibx-1];
00184             t2 = (1.0 - dely)*bvals[iby-1][ibxp1-1] + dely*bvals[ibyp1-1][ibxp1-1];
00185             dsky = avsky - (1.0 - delx)*t1 - delx*t2;
00186             map[kk+j] += dsky;
00187         }
00188     }
00189 
00190     /* Free some workspace */
00191 
00192     tidy();
00193     return(VIR_OK);
00194 }
00195 
00196 extern int imcore_backstats(ap_t *ap, float nullval, int satonly, 
00197                             float *skymed, float *skysig, float *sat) {
00198     int ilev,iclip,iloop,i,*ihist,isat,iter,*conf;
00199     long mpix,npix,k,mcpix,irej,lpix,nx,ny;
00200     float sigsq,skymedc,sigmac,*map,sata,fac;
00201 
00202     /* Get some info from the ap structure */
00203 
00204     map = ap->indata;
00205     conf = ap->confdata;
00206     nx = ap->lsiz;
00207     ny = ap->csiz;
00208 
00209     /* Check to make sure there are some non-zero values here */
00210 
00211     ilev = 1;
00212     for (i = 0; i < nx*ny; i++) {
00213         if (map[i] != 0.0 && conf[i] > 0) {
00214             ilev = 0;
00215             break;
00216         }
00217     }
00218     if (ilev == 1) {
00219         *skymed = 0.0;
00220         *skysig = 0.0;
00221         *sat = 0.0;
00222         return(VIR_WARN);
00223     }
00224 
00225     /* First, get some workspace for the background histogram */
00226 
00227     ihist = cpl_calloc(MAXHIST,sizeof(*ihist));
00228 
00229     /* Loop for up to 10 iterations. For each iteration we multiply the 
00230        input data by a successively higher power of 1 in order to 
00231        try and deal with data that has very small noise estimates */
00232 
00233     fac = 0.5;
00234     for (iter = 0; iter <=9; iter++) {
00235         fac *= 2.0;
00236         for (k = 0; k < MAXHIST; k++)
00237             ihist[k] = 0;
00238 
00239         /* Now form the histogram of all pixel intensities */
00240 
00241         mpix = 0;
00242         isat = 0;
00243         npix = nx*ny;
00244         for (k = 0; k < npix; k++) {
00245             if (map[k] != nullval && conf[k] > 0) {
00246                 ilev = MIN(MAXHISTVAL,MAX(MINHISTVAL,NINT(fac*map[k])));
00247                 ihist[ilev - MINHISTVAL] += 1;
00248                 isat = MAX(isat,ilev);
00249                 mpix++;
00250             }
00251         }
00252         sata = MIN(MAXHISTVAL,MAX(MINSATURATE,0.9*((float)isat))/fac);
00253         lpix = ihist[isat - MINHISTVAL];
00254         while (lpix < mpix/1000 && isat > MINHISTVAL) {
00255             isat--;
00256             lpix += ihist[isat - MINHISTVAL];
00257         }
00258         *sat = ((float)isat)/fac;
00259         *sat = MIN(MAXHISTVAL,MAX(MINSATURATE,MAX(0.95*(*sat),sata)));
00260 
00261         /* If all you want is the saturation level, then get out of here...*/
00262 
00263         if (satonly) {
00264             freespace(ihist);
00265             return(VIR_OK);
00266         }
00267 
00268         /* Now find the median and sigma */
00269 
00270         imcore_medsig(ihist,MAXHIST,MINHISTVAL-1,mpix,skymed,skysig);
00271         sigsq = (*skysig)*(*skysig);
00272 
00273         /* Do an iterative 3-sigma upper clip to give a more robust 
00274            estimator */
00275 
00276         iclip = MAXHISTVAL;
00277         mcpix = mpix;
00278         skymedc = *skymed;
00279         sigmac = *skysig;
00280         for (iloop = 0; iloop < 3; iloop++) {
00281             irej = 0;
00282             for (i = NINT(skymedc+3.0*sigmac); i <= iclip; i++)
00283                 irej += ihist[i - MINHISTVAL];
00284             if (irej == 0)
00285                 break;
00286             iclip = NINT(skymedc+3.0*sigmac)-1;
00287             mcpix = mcpix-irej;
00288             imcore_medsig(ihist,MAXHIST,MINHISTVAL-1,mcpix,&skymedc,&sigmac);
00289         }
00290         if (sigmac > 2.5)
00291             break;
00292     }
00293 
00294     /* Set the final answer */
00295 
00296     *skymed = skymedc/fac;
00297     *skysig = sigmac/fac;
00298     freespace(ihist);
00299     return(VIR_OK);
00300 }   
00301 
00302 extern void imcore_backest(ap_t *ap, float x, float y, float *skylev, 
00303                            float *skyrms) {
00304     int i,j,nbx,nby,nbsize,nbsizo2,iby,ibyp1,ibx,ibxp1;
00305     float **bvals,fnbsize,dely,delx,t1,t2;
00306 
00307     /* Define some local variables */
00308 
00309     nbx = ap->backmap.nbx;
00310     nby = ap->backmap.nby;
00311     nbsize = ap->backmap.nbsize;
00312     bvals = ap->backmap.bvals;
00313 
00314     /* Get closest pixel to the input location */
00315 
00316     i = NINT(x);
00317     j = NINT(y);
00318 
00319     /* Now, work out where in the map to do the interpolation */
00320 
00321     nbsizo2 = nbsize/2;
00322     fnbsize = 1.0/((float)nbsize);
00323     iby = (j + nbsizo2)/nbsize;
00324     ibyp1 = iby + 1;
00325     iby = MIN(nby,MAX(1,iby));
00326     ibyp1 = MIN(nby,ibyp1);
00327     dely = (j  - nbsize*iby + nbsizo2)*fnbsize;
00328     ibx = (i + nbsizo2)/nbsize;
00329     ibxp1 = ibx + 1;
00330     ibx = MIN(nbx,MAX(1,ibx));
00331     ibxp1 = MIN(nbx,ibxp1);
00332     delx = (i - nbsize*ibx + nbsizo2)*fnbsize;
00333 
00334     /* Now do a linear interpolation to find the background. Calculate MAD of
00335        the four adjacent background cells as an estimate of the RMS */
00336 
00337     t1 = (1.0 - dely)*bvals[iby-1][ibx-1] + dely*bvals[ibyp1-1][ibx-1];
00338     t2 = (1.0 - dely)*bvals[iby-1][ibxp1-1] + dely*bvals[ibyp1-1][ibxp1-1];
00339     *skylev = (1.0 - delx)*t1 + delx*t2;
00340     *skyrms = 0.25*(fabsf(bvals[iby-1][ibx-1] - *skylev) +
00341                     fabsf(bvals[ibyp1-1][ibx-1] - *skylev) +
00342                     fabsf(bvals[iby-1][ibxp1-1] - *skylev) +
00343                     fabsf(bvals[ibyp1-1][ibxp1-1] - *skylev));
00344 }
00345 
00346 extern void imcore_medsig(int *shist, int nh, int ist, int itarg, float *med, 
00347                           float *sig) {
00348     int isum, medata;
00349     float ffrac,sigmed;
00350  
00351     /* median */ 
00352 
00353     isum = 0;
00354     medata = ist;
00355     while (isum <= (itarg+1)/2 && (medata-MINHISTVAL) < nh) {
00356         medata++;
00357         isum += shist[medata-MINHISTVAL];
00358     }
00359     if (shist[medata-MINHISTVAL] == 0) {
00360         ffrac = 0.0;
00361     } else {
00362         ffrac = (float)(isum - (itarg+1)/2)/(float)shist[medata-MINHISTVAL];
00363     }
00364     *med = (float)medata - ffrac + 0.5;
00365 
00366     /* sigma */
00367 
00368     isum = 0;
00369     medata = ist;
00370     while (isum <= (itarg+3)/4 && (medata-MINHISTVAL) < nh) {
00371         medata++;
00372         isum += shist[medata-MINHISTVAL];
00373     }
00374     if (shist[medata-MINHISTVAL] == 0) {
00375         ffrac = 0.0;
00376     } else {
00377         ffrac = (float)(isum - (itarg+3)/4)/(float)shist[medata-MINHISTVAL];
00378     }
00379     sigmed = (float)medata - ffrac + 0.5;
00380     *sig = 1.48*(*med - sigmed);
00381     *sig = MAX(0.5,*sig);
00382 }
00383 
00384 static void sortit (float ia[], int n) {
00385     int i, j, ii, jj, ifin;
00386     float it;
00387  
00388     jj = 4;
00389     while (jj < n) 
00390         jj = 2 * jj;
00391     jj = MIN(n,(3 * jj)/4 - 1);
00392     while (jj > 1) {
00393         jj = jj/2;
00394         ifin = n - jj;
00395         for (ii = 0; ii < ifin; ii++) {
00396             i = ii;
00397             j = i + jj;
00398             if (ia[i] <= ia[j]) 
00399                 continue;
00400             it = ia[j];
00401             do {
00402                 ia[j] = ia[i];
00403                 j = i;
00404                 i = i - jj;
00405                 if (i < 0) 
00406                     break;
00407             } while (ia[i] > it);
00408             ia[j] = it;
00409         }
00410     }
00411     return;
00412 }
00413 
00414 
00415 static void tidy(void) {
00416     int i;
00417  
00418     freespace(nnp);
00419     if (hist != NULL) {
00420         for (i = 0; i < npvx; i++) 
00421             freespace(hist[i]);
00422     }
00423     freespace(hist);
00424     return;
00425 }
00426 
00427 /*
00428 
00429 $Log: imcore_background.c,v $
00430 Revision 1.8  2007/12/19 13:16:50  jim
00431 Modified the interation step size
00432 
00433 Revision 1.7  2007/04/23 12:51:27  jim
00434 Traps for 0 in confidence map
00435 
00436 Revision 1.6  2007/03/01 12:38:26  jim
00437 Small modifications after a bit of code checking
00438 
00439 Revision 1.5  2006/08/01 11:27:54  jim
00440 Modifications to imcore background estimation and to add ability to
00441 specify the smoothing kernel width
00442 
00443 Revision 1.4  2006/07/11 14:50:27  jim
00444 Modified to do arithmetic in double precision for stats
00445 
00446 Revision 1.3  2006/06/30 21:31:09  jim
00447 MOdifications to background routines and smoothing kernel
00448 
00449 Revision 1.2  2006/04/20 11:15:19  jim
00450 A few minor bugs fixed
00451 
00452 Revision 1.1  2005/09/13 13:25:28  jim
00453 Initial entry after modifications to make cpl compliant
00454 
00455 
00456 */
00457 

Generated on Sat Apr 6 04:03:07 2013 for VIRCAM Pipeline by  doxygen 1.5.1