imcore_conf.c

00001 /*
00002 
00003 $Id: imcore_conf.c,v 1.17 2007/05/03 11:15:34 jim Exp $
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     /* Initialise output */
00056 
00057     *outcat = NULL;
00058 
00059     /* Useful constants */
00060 
00061     fconst = 1.0/M_LN2;
00062     nullval = 0.0;
00063     cattype = cattyp;
00064     nobjects = 0;
00065 
00066     /* Open input image */
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     /* Open the associated confidence map, if it exists */
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     /* Get mflag array for flagging saturated pixels */
00097 
00098     npix = nx*ny;
00099     mflag = cpl_calloc(npix,sizeof(*mflag));
00100 
00101     /* Open the ap structure and define some stuff in it */
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     /* Open the output catalogue FITS table */
00120 
00121     tabinit(&ap);
00122 
00123     /* Compute a saturation level before background correction. */
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     /* Flag up regions where the value is above the saturation level
00130        determined above. */
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     /* Compute the background variation and remove it from the data*/
00139 
00140     retval = imcore_background(&ap,nbsize,nullval);
00141     if (retval != VIR_OK) 
00142         FATAL_ERR("Error calculating background");
00143 
00144     /* Compute a saturation level after background correction. */
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     /* Compute background statistics */
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     /* Get the propertly list for the input file and add some info*/
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     /* Take mean sky level out of data */
00168 
00169     for (i = 0; i < nx*ny; i++)
00170         indata[i] -= skymed;
00171     
00172     /* Work out isophotal detection threshold levels */
00173 
00174     thresh = threshold*skysig;
00175     
00176     /* Minimum intensity for consideration */
00177 
00178     xintmin = 1.5*thresh*((float)ipix);
00179 
00180     /* Minimum size for considering multiple images */
00181 
00182     mulpix = MAX(8,2*ipix);
00183 
00184     /* Actual areal profile levels: T, 2T, 4T, 8T,...but written wrt T
00185        i.e. threshold as a power of 2 */
00186 
00187     offset = logf(thresh)*fconst;
00188 
00189     /* Get a bit of workspace for buffers */
00190 
00191     smoothed = cpl_malloc(nx*sizeof(*smoothed));
00192     smoothedc = cpl_malloc(nx*sizeof(*smoothedc));
00193 
00194     /* Define a few things more things in ap structure */
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     /* Now update mflag array */
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     /* Set the weights */
00213 
00214     crweights(filtfwhm);
00215     nw2 = NW/2;
00216 
00217     /* Right, now for the extraction loop.  Begin by defining a group of
00218        three rows of data and confidence */
00219 
00220     for (j = nw2; j < ny-nw2; j++) {
00221         current = indata + j*nx;
00222         currentc = confdata + j*nx;
00223         convolve(j);
00224    
00225         /* Do the detection now */
00226 
00227         apline(&ap,current,currentc,smoothed,smoothedc,j,NULL);
00228 
00229         /* Make sure we're not overruning the stacks */
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         /* See if there are any images to terminate */
00238 
00239         if (ap.ipstack > 1)
00240             terminate(&ap);
00241     }
00242 
00243     /* Post process. First truncate the cpl_table to the correct size and then
00244        work out an estimate of the seeing */
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     /* Create a property list with extra parameters. First QC parameters */
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     /* Now DRS parameters */
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     /* Now wrap all this stuff up and send it back */
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     /* Tidy and exit */  
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     /* Get the kernel size */
00309 
00310     nw2 = NW/2;
00311     
00312     /* Set the normalisation constants */
00313 
00314     gsigsq = 1.0/(2.0*pow(MAX(1.0,(double)filtfwhm)/2.35,2.0));
00315     renorm = 0.0;
00316 
00317     /* Now work out the weights */
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     /* Now normalise the weights */
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     /* Zero the summations */
00351 
00352     for (i = 0; i < nx; i++) {
00353         smoothed[i] = 0.0;
00354         smoothedc[i] = 0.0;
00355     }
00356 
00357     /* Now big is the smoothing kernel? */
00358 
00359     nw2 = NW/2;
00360 
00361     /* Now loop for each column */
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 $Log: imcore_conf.c,v $
00390 Revision 1.17  2007/05/03 11:15:34  jim
00391 Fixed little problem with table wcs
00392 
00393 Revision 1.16  2007/05/02 09:11:35  jim
00394 Modified to allow for inclusion of table WCS keywords into FITS header
00395 
00396 Revision 1.15  2007/03/01 12:38:26  jim
00397 Small modifications after a bit of code checking
00398 
00399 Revision 1.14  2006/11/27 11:56:59  jim
00400 Changed cpl_propertylist_append to cpl_propertylist_update calls
00401 
00402 Revision 1.13  2006/08/01 11:27:54  jim
00403 Modifications to imcore background estimation and to add ability to
00404 specify the smoothing kernel width
00405 
00406 Revision 1.12  2006/07/11 14:51:02  jim
00407 Fixed small bug in the range of the main loop
00408 
00409 Revision 1.11  2006/07/03 09:33:19  jim
00410 Fixed a few things to keep the compiler happy
00411 
00412 Revision 1.10  2006/06/30 21:31:09  jim
00413 MOdifications to background routines and smoothing kernel
00414 
00415 Revision 1.9  2006/06/13 14:06:10  jim
00416 Gets gain estimate from header
00417 
00418 Revision 1.8  2006/06/06 13:04:22  jim
00419 Fixed apline so that it now takes confidence map info correctly
00420 
00421 Revision 1.7  2006/05/30 12:13:53  jim
00422 Initialised nobjects
00423 
00424 Revision 1.6  2006/05/18 12:35:01  jim
00425 Fixed bug for header writing
00426 
00427 Revision 1.5  2006/03/15 10:43:42  jim
00428 Fixed a few things
00429 
00430 Revision 1.4  2006/03/01 10:31:29  jim
00431 Now uses new vir_fits objects
00432 
00433 Revision 1.3  2005/09/22 08:40:42  jim
00434 Fixed some bugs in imcore and added classification Removed some unnecessary
00435 files
00436 
00437 Revision 1.2  2005/09/20 15:07:47  jim
00438 Fixed a few bugs and added a few things
00439 
00440 Revision 1.1  2005/09/13 13:25:29  jim
00441 Initial entry after modifications to make cpl compliant
00442 
00443 
00444 */

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