00001
00002
00003
00004
00005
00006
00007
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include <math.h>
00011
00012 #include <cpl.h>
00013 #include "../vircam_fits.h"
00014 #include "../vircam_pfits.h"
00015 #include "../vircam_wcsutils.h"
00016
00017 #include "ap.h"
00018 #include "util.h"
00019 #include "imcore.h"
00020 #include "floatmath.h"
00021
00022 #define FATAL_ERR(_a) {freetable(tab); cpl_msg_error(fctid,_a); tidy(); return(VIR_FATAL);}
00023
00024 #define NW 5
00025
00026 static float *smoothed = NULL;
00027 static float *smoothedc = NULL;
00028 static unsigned char *mflag = NULL;
00029 static float *indata = NULL;
00030 static int *confdata = NULL;
00031 static ap_t ap;
00032 static int freeconf = 0;
00033
00034 static float weights[NW*NW];
00035 static float weightc[NW*NW];
00036 static long nx;
00037 static long ny;
00038
00039 static void crweights(float);
00040 static void convolve(int);
00041 static void tidy(void);
00042
00043 extern int imcore_conf(vir_fits *infile, vir_fits *conf, int ipix,
00044 float threshold, int icrowd, float rcore, int nbsize,
00045 int cattyp, float filtfwhm, vir_tfits **outcat) {
00046
00047 int i,retval,mulpix,j,*currentc,nw2,status;
00048 float fconst,nullval,skymed,skysig,thresh,xintmin,offset;
00049 float *current,isatbc,isat,junk;
00050 long npix,nxc,nyc,npts;
00051 cpl_image *map,*cmap;
00052 cpl_propertylist *plist,*extra;
00053 const char *fctid = "imcore_conf";
00054
00055
00056
00057 *outcat = NULL;
00058
00059
00060
00061 fconst = 1.0/M_LN2;
00062 nullval = 0.0;
00063 cattype = cattyp;
00064 nobjects = 0;
00065
00066
00067
00068 map = vircam_fits_get_image(infile);
00069 if ((indata = cpl_image_get_data_float(map)) == NULL)
00070 FATAL_ERR("Error getting image data");
00071 nx = (long)cpl_image_get_size_x(map);
00072 ny = (long)cpl_image_get_size_y(map);
00073 npts = nx*ny;
00074 if (vircam_pfits_get_gain(vircam_fits_get_ehu(infile),&gain) != VIR_OK)
00075 gain = 5.0;
00076
00077
00078
00079 if (conf != NULL) {
00080 cmap = vircam_fits_get_image(conf);
00081 if ((confdata = cpl_image_get_data(cmap)) == NULL)
00082 FATAL_ERR("Error getting confidence map data");
00083 nxc = (long)cpl_image_get_size_x(cmap);
00084 nyc = (long)cpl_image_get_size_y(cmap);
00085 if ((nx != nxc) || (ny != nyc))
00086 FATAL_ERR("Input image and confidence dimensions don't match");
00087 freeconf = 0;
00088 } else {
00089 confdata = cpl_malloc(npts*sizeof(*confdata));
00090 for (i = 0; i < npts; i++)
00091 confdata[i] = 100;
00092 freeconf = 1;
00093 cmap = NULL;
00094 }
00095
00096
00097
00098 npix = nx*ny;
00099 mflag = cpl_calloc(npix,sizeof(*mflag));
00100
00101
00102
00103 ap.lsiz = nx;
00104 ap.csiz = ny;
00105 ap.inframe = map;
00106 ap.conframe = cmap;
00107 ap.xtnum = vircam_fits_get_nexten(infile);
00108 apinit(&ap);
00109 ap.indata = indata;
00110 ap.confdata = confdata;
00111 ap.multiply = 1;
00112 ap.ipnop = ipix;
00113 ap.mflag = mflag;
00114 ap.rcore = rcore;
00115 ap.filtfwhm = filtfwhm;
00116 ap.icrowd = icrowd;
00117 ap.fconst = fconst;
00118
00119
00120
00121 tabinit(&ap);
00122
00123
00124
00125 retval = imcore_backstats(&ap,nullval,1,&skymed,&skysig,&isatbc);
00126 if (retval != VIR_OK)
00127 FATAL_ERR("Error calculating saturation level");
00128
00129
00130
00131
00132 for (i = 0; i < npix ; i++)
00133 if (indata[i] > isatbc)
00134 mflag[i] = MF_SATURATED;
00135 else if (indata[i] < STUPID_VALUE)
00136 mflag[i] = MF_STUPID_VALUE;
00137
00138
00139
00140 retval = imcore_background(&ap,nbsize,nullval);
00141 if (retval != VIR_OK)
00142 FATAL_ERR("Error calculating background");
00143
00144
00145
00146 retval = imcore_backstats(&ap,nullval,1,&skymed,&skysig,&isat);
00147 if (retval != VIR_OK)
00148 FATAL_ERR("Error calculating saturation");
00149
00150
00151
00152
00153 retval = imcore_backstats(&ap,nullval,0,&skymed,&skysig,&junk);
00154 if (retval != VIR_OK)
00155 FATAL_ERR("Error calculating background stats");
00156
00157
00158
00159 plist = vircam_fits_get_ehu(infile);
00160 cpl_propertylist_update_float(plist,"ESO DRS SKYLEVEL",skymed);
00161 cpl_propertylist_set_comment(plist,"ESO DRS SKYLEVEL",
00162 "[adu] Median sky brightness");
00163 cpl_propertylist_update_float(plist,"ESO DRS SKYNOISE",skysig);
00164 cpl_propertylist_set_comment(plist,"ESO DRS SKYNOISE",
00165 "[adu] Pixel noise at sky level");
00166
00167
00168
00169 for (i = 0; i < nx*ny; i++)
00170 indata[i] -= skymed;
00171
00172
00173
00174 thresh = threshold*skysig;
00175
00176
00177
00178 xintmin = 1.5*thresh*((float)ipix);
00179
00180
00181
00182 mulpix = MAX(8,2*ipix);
00183
00184
00185
00186
00187 offset = logf(thresh)*fconst;
00188
00189
00190
00191 smoothed = cpl_malloc(nx*sizeof(*smoothed));
00192 smoothedc = cpl_malloc(nx*sizeof(*smoothedc));
00193
00194
00195
00196 ap.mulpix = mulpix;
00197 ap.areal_offset = offset;
00198 ap.thresh = thresh;
00199 ap.xintmin = xintmin;
00200 ap.sigma = skysig;
00201 ap.background = skymed;
00202 ap.saturation = (float)isat;
00203
00204
00205
00206 for (i = 0; i < npix; i++) {
00207 if ((indata[i]*confdata[i]*0.01 > 3.0*skysig) &&
00208 (mflag[i] != MF_SATURATED))
00209 mflag[i] = MF_UNSATPIX;
00210 }
00211
00212
00213
00214 crweights(filtfwhm);
00215 nw2 = NW/2;
00216
00217
00218
00219
00220 for (j = nw2; j < ny-nw2; j++) {
00221 current = indata + j*nx;
00222 currentc = confdata + j*nx;
00223 convolve(j);
00224
00225
00226
00227 apline(&ap,current,currentc,smoothed,smoothedc,j,NULL);
00228
00229
00230
00231 if (ap.ibstack > (ap.maxbl - ap.lsiz))
00232 apfu(&ap);
00233 if (ap.ipstack > (ap.maxpa*3/4))
00234 for (i = 0; i < ap.maxpa*3/8; i++)
00235 apfu(&ap);
00236
00237
00238
00239 if (ap.ipstack > 1)
00240 terminate(&ap);
00241 }
00242
00243
00244
00245
00246 cpl_table_set_size(tab,nobjects);
00247 retval = do_seeing(&ap);
00248 if (retval != VIR_OK)
00249 FATAL_ERR("Error working out seeing");
00250 tabclose(&ap);
00251
00252
00253
00254 extra = cpl_propertylist_duplicate(vircam_fits_get_ehu(infile));
00255 cpl_propertylist_update_float(extra,"ESO QC SATURATION",ap.saturation);
00256 cpl_propertylist_set_comment(extra,"ESO QC SATURATION",
00257 "[adu] Saturation level");
00258 cpl_propertylist_update_float(extra,"ESO QC MEAN_SKY",ap.background);
00259 cpl_propertylist_set_comment(extra,"ESO QC MEAN_SKY",
00260 "[adu] Median sky brightness");
00261 cpl_propertylist_update_float(extra,"ESO QC SKY_NOISE",ap.sigma);
00262 cpl_propertylist_set_comment(extra,"ESO QC SKY_NOISE",
00263 "[adu] Pixel noise at sky level");
00264
00265
00266
00267 cpl_propertylist_update_float(extra,"ESO DRS THRESHOL",ap.thresh);
00268 cpl_propertylist_set_comment(extra,"ESO DRS THRESHOL",
00269 "[adu] Isophotal analysis threshold");
00270 cpl_propertylist_update_int(extra,"ESO DRS MINPIX",ap.ipnop);
00271 cpl_propertylist_set_comment(extra,"ESO DRS MINPIX",
00272 "[pixels] Minimum size for images");
00273 cpl_propertylist_update_int(extra,"ESO DRS CROWDED",ap.icrowd);
00274 cpl_propertylist_set_comment(extra,"ESO DRS CROWDED",
00275 "Crowded field analysis flag");
00276 cpl_propertylist_update_float(extra,"ESO DRS RCORE",ap.rcore);
00277 cpl_propertylist_set_comment(extra,"ESO DRS RCORE",
00278 "[pixels] Core radius for default profile fit");
00279 cpl_propertylist_update_float(extra,"ESO DRS SEEING",ap.fwhm);
00280 cpl_propertylist_set_comment(extra,"ESO DRS SEEING",
00281 "[pixels] Average FWHM");
00282 cpl_propertylist_update_float(extra,"ESO DRS FILTFWHM",ap.filtfwhm);
00283 cpl_propertylist_set_comment(extra,"ESO DRS FILTFWHM",
00284 "[pixels] FWHM of smoothing kernel");
00285 cpl_propertylist_update_int(extra,"ESO DRS XCOL",imcore_xcol);
00286 cpl_propertylist_set_comment(extra,"ESO DRS XCOL","Column for X position");
00287 cpl_propertylist_update_int(extra,"ESO DRS YCOL",imcore_ycol);
00288 cpl_propertylist_set_comment(extra,"ESO DRS YCOL","Column for Y position");
00289
00290
00291
00292 plist = cpl_propertylist_duplicate(vircam_fits_get_phu(infile));
00293 status = VIR_OK;
00294 (void)vircam_tabwcs(extra,imcore_xcol,imcore_ycol,&status);
00295 *outcat = vircam_tfits_wrap(tab,NULL,plist,extra);
00296
00297
00298
00299 tidy();
00300 return(VIR_OK);
00301 }
00302
00303 static void crweights(float filtfwhm) {
00304 int i,j,nw2,n;
00305 double gsigsq,di,dj;
00306 float renorm;
00307
00308
00309
00310 nw2 = NW/2;
00311
00312
00313
00314 gsigsq = 1.0/(2.0*pow(MAX(1.0,(double)filtfwhm)/2.35,2.0));
00315 renorm = 0.0;
00316
00317
00318
00319 n = -1;
00320 for (i = -nw2; i <= nw2; i++) {
00321 di = (double)i;
00322 di *= gsigsq*di;
00323 for (j = -nw2; j <= nw2; j++) {
00324 dj = (double)j;
00325 dj *= gsigsq*dj;
00326 n++;
00327 weights[n] = (float)exp(-(di + dj));
00328 renorm += weights[n];
00329 }
00330 }
00331
00332
00333
00334 n = -1;
00335 for (i = -nw2; i <= nw2; i++) {
00336 for (j = -nw2; j <= nw2; j++) {
00337 n++;
00338 weights[n] /= renorm;
00339 weightc[n] = 0.01*weights[n];
00340 }
00341 }
00342 }
00343
00344 static void convolve(int ir) {
00345 int i,nw2,ix,jx,jy,n;
00346 float *idata;
00347 int *cdata;
00348
00349
00350
00351
00352 for (i = 0; i < nx; i++) {
00353 smoothed[i] = 0.0;
00354 smoothedc[i] = 0.0;
00355 }
00356
00357
00358
00359 nw2 = NW/2;
00360
00361
00362
00363 for (ix = nw2; ix < nx-nw2; ix++) {
00364 n = -1;
00365 for (jy = ir-nw2; jy <= ir+nw2; jy++) {
00366 idata = indata + jy*nx;
00367 cdata = confdata + jy*nx;
00368 for (jx = ix-nw2; jx <= ix+nw2; jx++) {
00369 n++;
00370 smoothed[ix] += weights[n]*idata[jx];
00371 smoothedc[ix] += weightc[n]*idata[jx]*cdata[jx];
00372 }
00373 }
00374 }
00375 }
00376
00377 static void tidy(void) {
00378
00379 if (freeconf)
00380 freespace(confdata);
00381 freespace(smoothed);
00382 freespace(smoothedc);
00383 freespace(mflag);
00384 apclose(&ap);
00385 }
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444