00001
00002
00003
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
00035
00036 map = ap->indata;
00037 conf = ap->confdata;
00038 mflag = ap->mflag;
00039 nx = ap->lsiz;
00040 ny = ap->csiz;
00041
00042
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));
00052
00053
00054
00055 nbx = nx/nbsize;
00056 nby = ny/nbsize;
00057 npixstripe = nbsize*nx;
00058 npvx = nbx;
00059 npvy = nby;
00060
00061
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
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
00074
00075 ap->backmap.nbx = nbx;
00076 ap->backmap.nby = nby;
00077 ap->backmap.nbsize = nbsize;
00078 ap->backmap.bvals = bvals;
00079
00080
00081
00082 nnp = cpl_malloc(nbx*sizeof(int));
00083
00084
00085
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
00093
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
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
00115
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
00140
00141 bfilt(bvals,nbx,nby);
00142
00143
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
00155
00156 nbsizo2 = nbsize/2;
00157 fnbsize = 1.0/((float)nbsize);
00158 for (k = 0; k < ny; k++) {
00159 kk = k*nx;
00160
00161
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
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
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
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
00203
00204 map = ap->indata;
00205 conf = ap->confdata;
00206 nx = ap->lsiz;
00207 ny = ap->csiz;
00208
00209
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
00226
00227 ihist = cpl_calloc(MAXHIST,sizeof(*ihist));
00228
00229
00230
00231
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
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
00262
00263 if (satonly) {
00264 freespace(ihist);
00265 return(VIR_OK);
00266 }
00267
00268
00269
00270 imcore_medsig(ihist,MAXHIST,MINHISTVAL-1,mpix,skymed,skysig);
00271 sigsq = (*skysig)*(*skysig);
00272
00273
00274
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
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
00308
00309 nbx = ap->backmap.nbx;
00310 nby = ap->backmap.nby;
00311 nbsize = ap->backmap.nbsize;
00312 bvals = ap->backmap.bvals;
00313
00314
00315
00316 i = NINT(x);
00317 j = NINT(y);
00318
00319
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
00335
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
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
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
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457