classify.c

00001 /*
00002 
00003 $Id: classify.c,v 1.9 2007/05/02 09:11:34 jim Exp $
00004 
00005 */
00006 
00007 #include <stdio.h>
00008 #include <stdlib.h>
00009 #include <math.h>
00010 #include <strings.h>
00011 
00012 #include <cpl.h>
00013 #include <vircam_utils.h>
00014 #include <vircam_pfits.h>
00015 #include "classify.h"
00016 
00017 #define MAXHIST   111
00018 #define STEP      0.05
00019 #define NSAMPLE   150
00020 #define MAXLOOP   5
00021 #define BLIMDEF   15.0;
00022 #define FLIMDEF   11.0;
00023 #define CMINDEF   7.5
00024 #define CMAXDEF   15.0
00025 #define NAREAL    8
00026 #define PI        3.14159265358979323846
00027 
00028 #define COREMAG(A,B,C) 2.5*log10((double)(max(C,A-B)))
00029 
00030 /* Make the data arrays and header values global */
00031 
00032 static long nrows;
00033 static float thresh,skylevel,skynoise,rcore,exptime;
00034 
00035 /* Derived values */
00036 
00037 static int poor;
00038 static float sigell,fitell,elllim,sigellf,fitellf;
00039 static float blim,flim,cmin,cmax;
00040 static float fit1,fit2,fit3,fit4,fit5;
00041 static float fit_final,sigma_final;
00042 static float *lower1,*lower2,*lower3,*upper1,*upper2,*upper3,*uppere;
00043 static float avsig1,avsig2,avsig3,wt1,wt2,wt3;
00044 
00045 /* Classification values */
00046 
00047 static int nstar,ngal,njunk,ncmp;
00048 
00049 /* Values for the data quality and aperture corrections */
00050 
00051 static float avsat,corlim,cormin,apcpkht,apcor,apcor1,apcor2,apcor3,apcor4;
00052 static float apcor5,apcor6,apcor7;
00053 
00054 /* Data arrays */
00055 
00056 static float *workspace = NULL;
00057 static cpl_table *catcopy = NULL;
00058 
00059 static float *core_flux,*core1_flux,*core2_flux,*core3_flux,*core4_flux;
00060 static float *core5_flux,*peak_height,*peak_mag,*ellipticity,*iso_flux;
00061 static float *total_flux,*cls,*sig;
00062 
00063 
00064 /* Column definitions */
00065 
00066 #define NCOL 10
00067 static const char *cols32[NCOL] = {"Core_flux","Core1_flux","Core2_flux",
00068                                    "Core3_flux","Core4_flux","Peak_height",
00069                                    "Ellipticity","Isophotal_flux","Total_flux",
00070                                  "Core5_flux"};
00071 
00072 static const char *cols80[NCOL] = {"Aper_flux_3","Aper_flux_1","Aper_flux_4",
00073                                    "Aper_flux_5","Aper_flux_6","Peak_height",
00074                                    "Ellipticity","Isophotal_flux","Hall_flux",
00075                                    "Aper_flux_7"};
00076 
00077 static int ncols;
00078 
00079 
00080 /* Subroutine prototypes */
00081 
00082 static void anhist(float *, int, float *, float *);
00083 static void boundries(float *, float *, float *, float, float, float, float,
00084                       int, float, float, float *, float *, float *, float *);
00085 static void boundpk(float *, float *, float, float, float *, float *,
00086                     float *, float *);
00087 static void classify_run(void);
00088 static void classstats(float *, float *, int, float *, float *);
00089 static void classstats_ap(float *, float *);
00090 static void classstats_el(void);
00091 static void classstats_ellf(float);
00092 static void classstats_final(void);
00093 static void medstat(float *, int, float *, float *);
00094 static void sort1(float *, int);
00095 static void sort2(float *, float *, int);
00096 
00097 extern int classify(vir_tfits *catalogue, cpl_propertylist *plist ,
00098                     int cattype) {
00099     float fwhm,*work,*areal[NAREAL];
00100     float pkht,ell,core,ap,delap,area,junk,arg;
00101     char *cols[NCOL],colname[32];
00102     cpl_propertylist *extra;
00103     cpl_table *cat;
00104     const char *fctid = "vircam_classify";
00105     int i,n,iap,i1,i2;
00106 
00107     /* Get some DQC info from the extra propertylist generated by imcore */
00108 
00109     extra = vircam_tfits_get_ehu(catalogue);
00110     thresh = cpl_propertylist_get_float(extra,"ESO DRS THRESHOL");
00111     skylevel = cpl_propertylist_get_float(extra,"ESO QC MEAN_SKY");
00112     skynoise = cpl_propertylist_get_float(extra,"ESO QC SKY_NOISE");
00113     rcore = cpl_propertylist_get_float(extra,"ESO DRS RCORE");
00114     fwhm = cpl_propertylist_get_float(extra,"ESO DRS SEEING");
00115 
00116     /* Get the exposure time */
00117 
00118     if (vircam_pfits_get_exptime(plist,&exptime) != VIR_OK) {
00119         cpl_msg_warning(fctid,"Unable to get expsoure time!");
00120         exptime = 10.0;
00121     }
00122     
00123     /* Get the number of columns and decide which column labels to use */
00124 
00125     cat = vircam_tfits_get_table(catalogue);
00126     ncols = cpl_table_get_ncol(cat);
00127     switch (cattype) {
00128     case 1:
00129         for (i = 0; i < NCOL; i++)
00130             cols[i] = (char *)cols32[i];
00131         break;
00132     case 2:
00133         for (i = 0; i < NCOL; i++)
00134             cols[i] = (char *)cols80[i];
00135         break;
00136     default:
00137         cpl_msg_error(fctid,"Don't recognise catalogues with %d columns: cattype == %d",ncols,cattype);
00138         return(VIR_FATAL);
00139     }
00140 
00141     /* Make a copy of the table as you are going to muck about with the
00142        column values. Get the column data */
00143 
00144     catcopy = cpl_table_duplicate(cat);
00145     nrows = cpl_table_get_nrow(cat);
00146     core_flux = cpl_table_get_data_float(catcopy,cols[0]);
00147     core1_flux = cpl_table_get_data_float(catcopy,cols[1]);
00148     core2_flux = cpl_table_get_data_float(catcopy,cols[2]);
00149     core3_flux = cpl_table_get_data_float(catcopy,cols[3]);
00150     core4_flux = cpl_table_get_data_float(catcopy,cols[4]);
00151     peak_height = cpl_table_get_data_float(catcopy,cols[5]);
00152     ellipticity = cpl_table_get_data_float(catcopy,cols[6]);
00153     iso_flux = cpl_table_get_data_float(catcopy,cols[7]);
00154     total_flux = cpl_table_get_data_float(catcopy,cols[8]);
00155     core5_flux = cpl_table_get_data_float(catcopy,cols[9]);
00156     cls = cpl_table_get_data_float(cat,"Classification");
00157     sig = cpl_table_get_data_float(cat,"Statistic");
00158 
00159     /* Get some workspace */
00160 
00161     workspace = cpl_malloc(2*nrows*sizeof(float));
00162     peak_mag = workspace;
00163     work = workspace + nrows;
00164     
00165     /* Convert fluxes to "magnitudes" */
00166 
00167     for (i = 0; i < nrows; i++) {
00168         core_flux[i] = COREMAG(core_flux[i],0.0,1.0);
00169         core1_flux[i] = COREMAG(core1_flux[i],0.0,1.0);
00170         core2_flux[i] = COREMAG(core2_flux[i],0.0,1.0);
00171         core3_flux[i] = COREMAG(core3_flux[i],0.0,1.0);
00172         core4_flux[i] = COREMAG(core4_flux[i],0.0,1.0);
00173         core5_flux[i] = COREMAG(core5_flux[i],0.0,1.0);
00174         iso_flux[i] = COREMAG(iso_flux[i],0.0,1.0);
00175         total_flux[i] = COREMAG(total_flux[i],0.0,1.0);
00176         peak_mag[i] = COREMAG(peak_height[i],skynoise,0.1);
00177     }
00178 
00179     /*  Now get the areal profile information. You'll need this in a sec */
00180 
00181     for (i = 0; i < NAREAL; i++) {
00182         sprintf(colname,"Areal_%d_profile",i+1);
00183         areal[i] = cpl_table_get_data_float(catcopy,colname);
00184     }
00185 
00186     /* What is the seeing like? */
00187 
00188     poor = 0;
00189     if (fwhm > max(5.0,rcore*sqrt(2.0)))
00190         poor = 1;
00191 
00192     /* Ok, now call the routine that does all the work */
00193 
00194     classify_run();
00195 
00196     /* Right, now get a better estimate of the seeing */
00197 
00198     n = 0;
00199     for (i = 0; i < nrows; i++) {
00200         pkht = peak_height[i];
00201         ell = ellipticity[i];
00202         core = core_flux[i];
00203         if (cls[i] == -1.0 && ell < elllim && core < corlim && 
00204             pkht > 10.0*thresh) { 
00205             ap = log(0.5*pkht/thresh)/log(2.0) + 1.0;
00206             iap = (int)ap;
00207             delap = ap - (float)iap;
00208             if (iap > 0 && iap < NAREAL && areal[1][i] > 0.0) {
00209                 i1 = (iap-1)*nrows + i;
00210                 i2 = iap*nrows + i;
00211                 area = areal[iap-1][i]*(1.0 - delap) + areal[iap][i]*delap;
00212                 work[n++] = 2.0*sqrt(area/PI);
00213             }
00214         }
00215     }
00216     if (n > 2) { 
00217         medstat(work,n,&fwhm,&junk);
00218        
00219         /* Allow for finite pixel size */
00220 
00221         arg = 0.25*PI*fwhm*fwhm - 1;
00222         fwhm = 2.0*sqrt(max(0.0,arg/PI));
00223        
00224     } else
00225         fwhm = -1.0;
00226 
00227     /* Tidy up a bit */
00228 
00229     freespace(workspace);
00230     freetable(catcopy);
00231 
00232     /* Write header results into extra property list. First the QC */
00233 
00234     cpl_propertylist_update_float(extra,"ESO QC IMAGE_SIZE",fwhm);
00235     cpl_propertylist_set_comment(extra,"ESO QC IMAGE_SIZE",
00236                                  "[pixels] Average FWHM of stellar objects");
00237     cpl_propertylist_update_float(extra,"ESO QC ELLIPTICITY",fitell);
00238     cpl_propertylist_set_comment(extra,"ESO QC ELLIPTICITY",
00239                                  "Average stellar ellipticity (1-b/a)");
00240     switch (cattype) {
00241     case 1:
00242         cpl_propertylist_update_float(extra,"ESO QC APERTURE_CORR",apcor);
00243         cpl_propertylist_set_comment(extra,"ESO QC APERTURE_CORR",
00244                                      "Stellar ap-corr 1x core flux");
00245         break;
00246     case 2:
00247         cpl_propertylist_update_float(extra,"ESO QC APERTURE_CORR",apcor3);
00248         cpl_propertylist_set_comment(extra,"ESO QC APERTURE_CORR",
00249                                      "Stellar ap-corr 1x core flux");
00250         break;
00251     }
00252     cpl_propertylist_update_int(extra,"ESO QC NOISE_OBJ",njunk);
00253     cpl_propertylist_set_comment(extra,"ESO QC NOISE_OBJ",
00254                                  "Number of noise objects");
00255     cpl_propertylist_update_float(extra,"ESO QC SATURATION",avsat);
00256     
00257     /* Now some helpful DRS keywords */
00258 
00259     cpl_propertylist_update_bool(extra,"ESO DRS CLASSIFD",1);
00260     cpl_propertylist_set_comment(extra,"ESO DRS CLASSIFD",
00261                                  "Catalogue has been classified");
00262 
00263     /* Now the aperture correction keywords */
00264 
00265     switch (cattype) {
00266     case 1:
00267         cpl_propertylist_update_float(extra,"APCORPK",apcpkht);
00268         cpl_propertylist_set_comment(extra,"APCORPK","Stellar aperture correction - peak height");
00269         cpl_propertylist_update_float(extra,"APCOR1",apcor1);
00270         cpl_propertylist_set_comment(extra,"APCOR1","Stellar aperture correction - 1/2x core flux");
00271         cpl_propertylist_update_float(extra,"APCOR",apcor);
00272         cpl_propertylist_set_comment(extra,"APCOR","Stellar aperture correction - 1x core flux");
00273         cpl_propertylist_update_float(extra,"APCOR2",apcor2);
00274         cpl_propertylist_set_comment(extra,"APCOR2","Stellar aperture correction - sqrt(2)x core flux");
00275         cpl_propertylist_update_float(extra,"APCOR3",apcor3);
00276         cpl_propertylist_set_comment(extra,"APCOR3","Stellar aperture correction - 2x core flux");
00277         cpl_propertylist_update_float(extra,"APCOR4",apcor4);
00278         cpl_propertylist_set_comment(extra,"APCOR4","Stellar aperture correction - 2*sqrt(2)x core flux");
00279         cpl_propertylist_update_float(extra,"APCOR5",apcor5);
00280         cpl_propertylist_set_comment(extra,"APCOR5","Stellar aperture correction - 4x core flux");
00281         break;
00282     case 2:
00283         cpl_propertylist_update_float(extra,"APCORPK",apcpkht);
00284         cpl_propertylist_set_comment(extra,"APCORPK","Stellar aperture correction - peak height");
00285         cpl_propertylist_update_float(extra,"APCOR1",apcor1);
00286         cpl_propertylist_set_comment(extra,"APCOR1","Stellar aperture correction - 1/2x core flux");
00287         cpl_propertylist_update_float(extra,"APCOR2",apcor2);
00288         cpl_propertylist_set_comment(extra,"APCOR2","Stellar aperture correction - core/sqrt(2) flux");
00289         cpl_propertylist_update_float(extra,"APCOR3",apcor3);
00290         cpl_propertylist_set_comment(extra,"APCOR3","Stellar aperture correction - 1x core flux");
00291         cpl_propertylist_update_float(extra,"APCOR4",apcor4);
00292         cpl_propertylist_set_comment(extra,"APCOR4","Stellar aperture correction - sqrt(2)x core flux");
00293         cpl_propertylist_update_float(extra,"APCOR5",apcor5);
00294         cpl_propertylist_set_comment(extra,"APCOR5","Stellar aperture correction - 2x core flux");
00295         cpl_propertylist_update_float(extra,"APCOR6",apcor6);
00296         cpl_propertylist_set_comment(extra,"APCOR6","Stellar aperture correction - 2*sqrt(2)x core flux");
00297         cpl_propertylist_update_float(extra,"APCOR7",apcor7);
00298         cpl_propertylist_set_comment(extra,"APCOR7","Stellar aperture correction - 4x core flux");
00299 
00300         break;
00301     }
00302 
00303     /* Write header information to help GAIA */
00304 
00305     cpl_propertylist_update_string(extra,"SYMBOL1","{Ellipticity Position_angle Areal_1_profile Classification} {el");
00306     cpl_propertylist_update_string(extra,"SYMBOL2","lipse blue (1.0-$Ellipticity) $Position_angle+90 {} $Classific");
00307     cpl_propertylist_update_string(extra,"SYMBOL3","ation==1} {sqrt($Areal_1_profile*(1.0-$Ellipticity)/3.142)} : {");
00308     cpl_propertylist_update_string(extra,"SYMBOL4","Ellipticity Position_angle Areal_1_profile Classification} {el");
00309     cpl_propertylist_update_string(extra,"SYMBOL5","lipse red (1.0-$Ellipticity) $Position_angle+90 {} $Classific");
00310     cpl_propertylist_update_string(extra,"SYMBOL6","ation==-1} {sqrt($Areal_1_profile*(1.0-$Ellipticity)/3.142)} :");
00311     cpl_propertylist_update_string(extra,"SYMBOL7","{Ellipticity Position_angle Areal_1_profile Classification} {el");
00312     cpl_propertylist_update_string(extra,"SYMBOL8","lipse green (1.0-$Ellipticity) $Position_angle+90 {} $Classifi");
00313     cpl_propertylist_update_string(extra,"SYMBOL9","cation==0} {sqrt($Areal_1_profile*(1.0-$Ellipticity)/3.142)}");
00314 
00315     /* Get out of here */
00316 
00317     return(VIR_OK);
00318 }
00319 
00320 
00321 static void anhist(float *data, int n, float *medval, float *sigma) {
00322     int i,*histo,ilev,imax,ismax;
00323     float *sval,hmax,smax,hlim,ratio;
00324 
00325     /* Get some workspace for the histogram */
00326 
00327     histo = cpl_calloc(MAXHIST,sizeof(int));
00328     sval = cpl_calloc(MAXHIST,sizeof(float));
00329 
00330     /* Sort data into the histogram */
00331 
00332     for (i = 0; i < n; i++) {
00333         ilev = vircam_nint(data[i]/STEP);
00334         if (ilev >= -10 && ilev <= 100) {
00335             ilev += 10;
00336             histo[ilev] += 1;
00337         }
00338     }
00339 
00340     /* Now find the maximum of the histogram and its position... */
00341 
00342     hmax = 0.0;
00343     imax = 0;
00344     for (i = 0; i < MAXHIST; i++) {
00345         if (histo[i] > hmax) {
00346             hmax = (float)histo[i];
00347             imax = i;
00348         }
00349     }
00350 
00351     /* Trap for hmax == 0 */
00352 
00353     if (hmax == 0.0) {
00354         if (n >= 10) {
00355             *medval = data[(n+1)/2-1];
00356             *sigma = 1.48*0.5*(data[(3*n+3)/4-1] - data[(n+3)/4-1]);
00357         } else {
00358             *medval = 0.0;
00359             *sigma = 1.0;
00360         }
00361         return;
00362     }
00363 
00364     /* Now do three point running average to see if there are other local
00365        maxima */
00366 
00367     smax = 0.0;
00368     ismax = 0;
00369     for (i = 1; i < MAXHIST-1; i++) {
00370         sval[i] = (histo[i-1] + histo[i] + histo[i+1])/3.0;
00371         if (sval[i] > smax) {
00372             smax = sval[i];
00373             ismax = i;
00374         }
00375     }
00376     if (ismax < imax) {
00377         imax = ismax;
00378         hmax = (float)histo[imax];
00379     }
00380 
00381     /* Now check for lower local maxima */
00382 
00383     for (i = imax-1; i > 0; i--) {
00384         if (sval[i] >= sval[i+1] && sval[i] >= sval[i-1]) {
00385             if (sval[i] > 0.5*smax)
00386                 ismax = i;
00387         }
00388     }
00389     if (ismax < imax) {
00390         imax = ismax;
00391         hmax = (float)histo[imax];
00392     }
00393     
00394     /* Now work out where the peak is */
00395     
00396     *medval = min((float)(imax-10)*STEP,data[(n+1)/2-1]);
00397     hlim = vircam_nint(0.5*hmax);
00398     i = 1;
00399     while (histo[imax-i] > hlim && imax-i > 1)
00400         i++;
00401     ratio = hmax/max(1.0,(float)histo[imax-i]);
00402     *sigma = (float)i*STEP/(sqrt(2.0)*max(1.0,log(ratio)));
00403     *sigma = max(*sigma,0.5*STEP);
00404 
00405     /* Tidy and exit */
00406     
00407     freespace(histo);
00408     freespace(sval);
00409 }
00410 
00411 
00412 static void boundries(float *core1, float *core2, float *core3, float medval1,
00413                       float sigma1, float medval2, float sigma2, int small, 
00414                       float area1, float area2, float *wt, float *avsig, 
00415                       float *lower, float *upper) {
00416     int i,n;
00417     float c1,c2,dc,*work,xnoise,xmag,xflux,ratio,asign,junk;
00418 
00419     /* Get a workspace */
00420 
00421     work = cpl_malloc(nrows*sizeof(float));
00422 
00423     /* Initialise the lower boundry */
00424 
00425     lower[0] = cmin;
00426     lower[1] = cmax;
00427     asign = ((small == 1) ? -1.0 : 1.0);
00428     
00429     /* Now collect the data */
00430 
00431     n = 0;
00432     for (i = 0; i < nrows; i++) {
00433         c1 = core1[i];
00434         if (! poor) {
00435             c2 = core2[i];
00436             dc = asign*(c2 - c1);
00437             if (dc > medval1 - 3.0*sigma1 && c1 < 12.0)
00438                 work[n++] = dc - medval1;
00439         } else {
00440             c2 = core3[i];
00441             dc = c2 - c1;
00442             if (dc > medval2 - 3.0*sigma2 && c1 < 12.0)
00443                 work[n++] = dc - medval2;
00444         }
00445     }
00446  
00447     /* Find the median */
00448 
00449     medstat(work,n,avsig,&junk);
00450     freespace(work);
00451 
00452     /* Work out sigma levels for both types of seeing */
00453 
00454     if (! poor) {
00455         *wt = min(5.0,max(1.0,*avsig/sigma1));
00456         xnoise = sqrt(area1)*skynoise;
00457     } else {
00458         *wt = min(2.5,max(1.0,*avsig/sigma2));
00459         xnoise = sqrt(area2)*skynoise;
00460     }
00461 
00462     /* Now work out the boundries */
00463 
00464     for (i = 0; i < NSAMPLE; i++) {
00465         xmag = 5.0 + (float)(i+1)*0.1;
00466         xflux = pow(10.0,(double)(0.4*xmag));
00467         ratio = COREMAG(1.0+xnoise/xflux,0.0,0.0);
00468         if (! poor) {
00469             lower[i] = medval1 - 3.0*sqrt(sigma1*sigma1 + ratio*ratio);
00470             upper[i] = medval1 + 3.0*sqrt(sigma1*sigma1 + 0.5*ratio*ratio);
00471         } else {
00472             lower[i] = medval2 - 3.0*sqrt(sigma2*sigma2 + ratio*ratio);
00473             upper[i] = medval2 + 3.0*sqrt(sigma2*sigma2 + 0.5*ratio*ratio);
00474         }
00475     }
00476     upper[0] = ((poor == 0) ? medval1 : medval2);
00477     upper[1] = upper[0];
00478 }
00479 
00480 static void boundpk(float *core, float *pkht, float medval, float sigma, 
00481                     float *wt, float *avsig, float *lower, float *upper) {
00482     int i,n;
00483     float c,p,*work,xnoise,xmag,pmag,xflux,pflux,ratio,junk;
00484 
00485     /* Get the space for the boundry lines and a workspace */
00486 
00487     work = cpl_malloc(nrows*sizeof(float));
00488 
00489     /* Collect the data */
00490 
00491     n = 0;
00492     for (i = 0; i < nrows; i++) {
00493         c = core[i];
00494         p = pkht[i];
00495         if (c - p > medval - 3.0*sigma && c < 12)
00496             work[n++] = c - p - medval;
00497     }
00498 
00499     /* Find the median */
00500 
00501     medstat(work,n,avsig,&junk);
00502     freespace(work);
00503     *wt = min(5.0,max(1.0,*avsig/sigma));
00504 
00505     /* Now work out boundries */
00506 
00507     xnoise = sqrt(PI*rcore*rcore)*skynoise;
00508     for (i = 0; i < NSAMPLE; i++) {
00509         xmag = 5.0 + (float)(i+1)*0.1;
00510         pmag = xmag - medval;
00511         xflux = pow(10.0,(double)(0.4*xmag));
00512         pflux = pow(10.0,(double)(0.4*pmag));
00513         ratio = 2.5*log10((double)(1.0+max(xnoise/xflux,skynoise/pflux)));
00514         lower[i] = medval - 3.0*sqrt(sigma*sigma + ratio*ratio);
00515         upper[i] = medval + 3.0*sqrt(sigma*sigma + 0.5*ratio*ratio);
00516     }
00517     upper[0] = medval;
00518     upper[1] = upper[0];
00519 }
00520 
00521 
00522 static void classify_run(void) {
00523     float fluxlim,ell,pk,pkht,core,sig1,sig2,sig3,denom,w1,w2,w3;
00524     float core_small,core_large,core_midd,statistic,statcut,sigtot;
00525     float fit0,sigma0,xnoise,xmag,ratio,xflux,ratell,ratscl,ellbound;
00526     float *lower,*upper,sigma1,sigma2,sigma3,sigma4,sigma5;
00527     int i,iarg;
00528 
00529     /* Update faint limit to cope with short exposures */
00530 
00531     blim = BLIMDEF;
00532     flim = FLIMDEF;
00533     fluxlim = 2.5*log10((double)(5.0*sqrt(PI*rcore*rcore)*skynoise));
00534     flim = min(flim,max(6.0,fluxlim+3.0));
00535     corlim = min(blim,max(12.5,fluxlim+7.0));
00536     cormin = min(blim,max(12.5,fluxlim+6.5));
00537 
00538     /* Work out min and max core flux */
00539 
00540     cmin = CMINDEF;
00541     cmax = CMAXDEF;
00542     for (i = 0; i < nrows; i++) {
00543         xflux = core_flux[i];
00544         cmin = min(cmin,xflux);
00545         cmax = max(cmax,xflux);
00546     }
00547     cmin = max(fluxlim-0.5,cmin);
00548     cmax += 0.1;
00549     cmax = min(cmax,20.0);
00550 
00551     /* Work out ellipticity stats for likely stellar objects */
00552 
00553     classstats_el();
00554 
00555     /* Ok, get the classification statistics for each of the tests.  First
00556        the core flux vs 1/2*core flux */
00557 
00558     classstats(core_flux,core1_flux,1,&fit1,&sigma1);
00559 
00560     /* Core flux vs 2*core flux */
00561 
00562     classstats(core_flux,core3_flux,0,&fit2,&sigma2);
00563 
00564     /* Core flux vs sqrt(2)*core flux */
00565 
00566     classstats(core_flux,core2_flux,0,&fit4,&sigma4);
00567 
00568     /* Core flux vs 2*sqrt(2)*core flux */
00569 
00570     classstats(core_flux,core4_flux,0,&fit5,&sigma5);
00571 
00572     /* Core flux vs Peak height */
00573 
00574     classstats(core_flux,peak_mag,1,&fit3,&sigma3);
00575 
00576     /* Faint end ellipticity */
00577 
00578     classstats_ellf(fluxlim);
00579 
00580     /* Get workspace for the boundry arrays */
00581    
00582     lower1 = cpl_malloc(NSAMPLE*sizeof(float));
00583     lower2 = cpl_malloc(NSAMPLE*sizeof(float));
00584     lower3 = cpl_malloc(NSAMPLE*sizeof(float));
00585     upper1 = cpl_malloc(NSAMPLE*sizeof(float));
00586     upper2 = cpl_malloc(NSAMPLE*sizeof(float));
00587     upper3 = cpl_malloc(NSAMPLE*sizeof(float));
00588 
00589     /* Right, work out the boundries for the classification tests 
00590        First core vs sqrt(2)*core or core vs 0.5*core depending upon
00591        the seeing */
00592 
00593     boundries(core_flux,core1_flux,core2_flux,fit1,sigma1,fit4,sigma4,
00594               1,PI*rcore*rcore,2.0*PI*rcore*rcore,&wt1,&avsig1,lower1,
00595               upper1);
00596 
00597     /* Now core vs 2*core or core vs 2*sqrt(2)*core */
00598 
00599     boundries(core_flux,core3_flux,core4_flux,fit2,sigma2,fit5,sigma5,
00600               0,4.0*PI*rcore*rcore,8.0*PI*rcore*rcore,&wt2,&avsig2,lower2,
00601               upper2);
00602 
00603     /* Now core vs peak height */
00604 
00605     boundpk(core_flux,peak_mag,fit3,sigma3,&wt3,&avsig3,lower3,upper3); 
00606 
00607      
00608     /* Do final classification statistics and find the saturation limit */
00609 
00610     classstats_final();
00611 
00612     /* Define final boundries */
00613 
00614     lower = cpl_malloc(NSAMPLE*sizeof(float));
00615     upper = cpl_malloc(NSAMPLE*sizeof(float));
00616     uppere = cpl_malloc(NSAMPLE*sizeof(float));
00617     xnoise = sqrt(PI*rcore*rcore)*skynoise;
00618     ratell = xnoise/pow(10.0,0.4*(fluxlim+1.5));
00619     ratell = COREMAG(1.0+ratell,0.0,0.0);
00620     ratscl = (pow((fitellf + 2.0*sigellf - fitell),2.0) - 4.0*sigell*sigell)/(4.0*ratell*ratell);
00621     ratscl = max(0.25,min(10.0,ratscl));
00622     for (i = 0; i < NSAMPLE; i++) {
00623         xmag = 5.0 + 0.1*(float)(i+1);
00624         xflux = pow(10.0,0.4*xmag);
00625         ratio = 2.5*log10(1.0+xnoise/xflux);
00626         lower[i] = fit_final - 5.0*sqrt(sigma_final*sigma_final + ratio*ratio);
00627         upper[i] = fit_final + sqrt(9.0*sigma_final*sigma_final + 0.0*ratio*ratio);
00628         uppere[i] = fitell + 2.0*sqrt(sigell*sigell + ratscl*ratio*ratio);
00629         uppere[i] = min(0.5,uppere[i]);
00630     }
00631     elllim = min(0.5,max(0.25,fitell+2.0*sigell));
00632     fluxlim = 2.5*log10((double)(2.5*sqrt(PI*rcore*rcore)*skynoise));
00633 
00634     /* Ok, final classification loop now... */
00635 
00636     nstar = 0;
00637     ngal = 0;
00638     njunk = 0;
00639     ncmp = 0;
00640     for (i = 0; i < nrows; i++) {
00641         ell = ellipticity[i];
00642         pk = peak_height[i] + skylevel;
00643         pkht = peak_mag[i];
00644         core = core_flux[i];
00645         iarg = vircam_nint(10.0*(core - 5.0));
00646         iarg = max(1,min(NSAMPLE,iarg)) - 1;
00647         if (! poor) {
00648             sig1 = max(0.01,(fit1 - lower1[iarg])/3.0);
00649             sig2 = max(0.01,(fit2 - lower2[iarg])/3.0);
00650         } else {
00651             sig1 = max(0.01,(fit4 - lower1[iarg])/3.0);
00652             sig2 = max(0.01,(fit5 - lower2[iarg])/3.0);
00653         }
00654         sig3 = max(0.01,(fit3 - lower3[iarg])/3.0);
00655         denom = (wt1/sig1 + wt2/sig2 + wt3/sig3);
00656         w1 = (wt1/sig1)/denom;
00657         w2 = (wt2/sig2)/denom;
00658         w3 = (wt3/sig3)/denom;
00659         if (! poor) {
00660             core_small = core1_flux[i];
00661             core_large = core3_flux[i];
00662             statistic = (core - core_small - fit1)*w1 + 
00663                 (max(-3.0*sig2,core_large - core - fit2))*w2 + 
00664                 (core - pkht - fit3)*w3;
00665         } else {
00666             core_midd = core2_flux[i];
00667             core_large = core4_flux[i];
00668             statistic = (core_midd - core - fit4)*w1 +
00669                 (max(-3.0*sig2,core_large - core - fit5))*w2 + 
00670                 (core - pkht - fit3)*w3;
00671         }
00672         cls[i] = -1.0;
00673         statcut = upper[iarg] + 3.0*sigma_final*(exp(max(0.0,core-corlim+1.0)) - 1.0);
00674         if (statistic  >= statcut) 
00675             cls[i] = 1.0;
00676         else if (statistic <= lower[iarg])
00677             cls[i] = 0.0;
00678 
00679         /* Save distance from the stellar locus */
00680 
00681         sigtot = (fit_final - lower[iarg])/5.0;
00682         sig[i] = (statistic - fit_final)/sigtot;
00683 
00684         /* Right, now here are lots of overrides for special circumstances */
00685         /* Too spikey? -> junk */
00686 
00687         if (core - pkht - fit3 < -4.0*sig3) 
00688             cls[i] = 0.0;
00689 
00690         /* Elliptical star? -> compact */
00691 
00692         ellbound = max(elllim,uppere[iarg]);
00693         if (ell > ellbound && cls[i] == -1.0 && core < flim && sig[i] > -2.0)
00694             cls[i] = -2.0;
00695 
00696         /* Saturated? -> star */
00697 
00698         if (core > corlim && statistic >= lower[iarg])
00699             cls[i] = -1.0;
00700 
00701         /* Too elliptical? -> junk */
00702 
00703         if (ell > 0.9 && core < corlim)
00704             cls[i] = 0.0;
00705 
00706         /* Too faint? -> junk */
00707 
00708         if (core < fluxlim)
00709             cls[i] = 0.0;
00710 
00711         /* Now count how many you have of each */
00712 
00713         if (cls[i] == -1.0)
00714             nstar++;
00715         else if (cls[i] == 1.0)
00716             ngal++;
00717         else if (cls[i] == -2.0)
00718             ncmp++;
00719         else
00720             njunk++;
00721     }
00722 
00723     /* Do stats to get the aperture corrections */
00724 
00725     classstats_ap(&fit0,&sigma0);
00726     fit0 = max(fit5,fit0);
00727     apcpkht = fit0 + fit3; /* pkht */
00728     switch (ncols) {
00729     case 32:
00730         apcor1 = fit0 + fit1;  /* 0.5*core */
00731         apcor = fit0;          /* core */
00732         apcor2 = fit0 - fit4;  /* sqrt(2)*core */
00733         apcor3 = fit0 - fit2;  /* 2*core */
00734         apcor4 = fit0 - fit5;  /* 2*sqrt(2)*core */
00735         apcor5 = 0.0;          /* 4*core */
00736         break;
00737     case 80:
00738         apcor1 = fit0 + fit1;      /* 0.5*core */
00739         apcor2 = fit0 + 0.5*fit1;  /* 1/sqrt(2) * core */
00740         apcor3 = fit0;             /* core */
00741         apcor4 = fit0 - fit4;      /* core * sqrt(2) */
00742         apcor5 = fit0 - fit2;      /* 2*core */
00743         apcor6 = fit0 - fit5;      /* 2*sqrt(2)*core */
00744         apcor7 = 0.0;              /* 4*core */
00745         break;
00746     }
00747 
00748     /* Ok, now get rid of some workspace */
00749         
00750     freespace(lower1);
00751     freespace(lower2);
00752     freespace(lower3);
00753     freespace(upper1);
00754     freespace(upper2);
00755     freespace(upper3);
00756     freespace(lower);
00757     freespace(upper);
00758     freespace(uppere);
00759 
00760 }
00761 
00762 static void classstats(float *core1, float *core2, int small, float *medval, 
00763                        float *sigma) {
00764 
00765     int i,iloop,n;
00766     float *work,*dc,sigmaold,amult;
00767 
00768     /* Initialise the output values to something stupid */
00769 
00770     *medval = 0.0;
00771     *sigma = 1.0e6;
00772     amult = (small == 1 ? -1.0 : 1.0);
00773 
00774     /* Get some workspace */
00775 
00776     work = cpl_malloc(nrows*sizeof(float));
00777     dc = cpl_malloc(nrows*sizeof(float));
00778 
00779     /* Work out differences */
00780 
00781     for (i = 0; i < nrows; i++)
00782         dc[i] = amult*(core2[i] - core1[i]);
00783 
00784     /* Do an iteration loop */
00785 
00786     for (iloop = 0; iloop < MAXLOOP; iloop++) {
00787         sigmaold = *sigma;
00788         n = 0;
00789 
00790         /* Ok, gather up all the stats */
00791 
00792         for (i = 0; i < nrows; i++) {
00793             
00794             /* Clipping criteria */
00795 
00796             if (ellipticity[i] < elllim && core1[i]< blim && core1[i] > flim &&
00797                 fabs(dc[i] - *medval) < 3.0*(*sigma)) {
00798                 work[n++] = dc[i];
00799             }
00800         }
00801 
00802         /* Sort the work array and find the median and sigma */
00803 
00804         if (n > 2) {
00805             sort1(work,n);
00806             if (iloop == 0 && n > 10) {
00807                 anhist(work,n,medval,sigma);
00808             } else {
00809                 medstat(work,n,medval,sigma);
00810                 *sigma = min(sigmaold,*sigma);
00811             }
00812         } else {
00813             *medval = 0.0;
00814             *sigma = 0.01;
00815         }
00816 
00817         /* Just in case... */
00818 
00819         *sigma = max(*sigma,0.01);
00820     }
00821 
00822     /* Tidy and exit */
00823 
00824     freespace(work);
00825     freespace(dc);
00826 }
00827 
00828 static void classstats_el(void) {
00829     int iloop,n,i;
00830     float *work;
00831 
00832     /* Initialise the mean and sigma to something stupid */
00833 
00834     sigell = 1.0e6;
00835     fitell = 0.0;
00836 
00837     /* Get some workspace */
00838 
00839     work = cpl_malloc(nrows*sizeof(float));
00840 
00841     /* Do iteration loop */
00842 
00843     for (iloop = 0; iloop < MAXLOOP; iloop++) {
00844         n = 0;
00845         for (i = 0; i < nrows; i++) {
00846             if (ellipticity[i] < 0.5 && core_flux[i] < blim && core_flux[i] > flim &&
00847                 fabs(ellipticity[i] - fitell) < 2.0*sigell)
00848                 work[n++] = ellipticity[i];
00849         }
00850         if (n > 2)
00851             medstat(work,n,&fitell,&sigell);
00852         else {
00853             fitell = 0.25;
00854             sigell = 0.05;
00855         }
00856     }
00857     elllim = min(0.5,max(0.2,fitell+2.0*sigell));
00858 
00859     /* Get out of here */
00860 
00861     freespace(work);
00862 }
00863 
00864 static void classstats_ellf(float fluxlim) {
00865     int iloop,n,i;
00866     float *work;
00867 
00868     /* Initialise the mean and sigma to something stupid */
00869 
00870     sigellf = 1.0e6;
00871     fitellf = 0.0;
00872 
00873     /* Get some workspace */
00874 
00875     work = cpl_malloc(nrows*sizeof(float));
00876 
00877     /* Do iteration loop */
00878 
00879     for (iloop = 0; iloop < MAXLOOP; iloop++) {
00880         n = 0;
00881         for (i = 0; i < nrows; i++) {
00882             if (ellipticity[i] < 0.75 && core_flux[i] > fluxlim+1.0 && 
00883                 core_flux[i] < fluxlim+2.0 &&
00884                 fabs(ellipticity[i] - fitellf) < 2.0*sigellf)
00885                 work[n++] = ellipticity[i];
00886         }
00887         if (n > 2)
00888             medstat(work,n,&fitellf,&sigellf);
00889         else {
00890             fitellf = 0.25;
00891             sigellf = 0.05;
00892         }
00893     }
00894 
00895     /* Get out of here */
00896 
00897     freespace(work);
00898 }
00899 
00900 static void classstats_ap(float *medval, float *sigma) {
00901 
00902     int i,iloop,n;
00903     float *work,*dc,c2;
00904 
00905     /* Initialise the output values to something stupid */
00906 
00907     *medval = 0.0;
00908     *sigma = 1.0e6;
00909     elllim = min(0.5,max(0.2,fitell+2.0*sigell));
00910 
00911     /* Get some workspace */
00912 
00913     work = cpl_malloc(nrows*sizeof(float));
00914     dc = cpl_malloc(nrows*sizeof(float));
00915 
00916     /* Work out differences */
00917 
00918     for (i = 0; i < nrows; i++) {
00919         c2 = max(iso_flux[i],max(total_flux[i],core5_flux[i]));
00920         dc[i] = c2 - core_flux[i];
00921     }
00922 
00923     /* Do an iteration loop */
00924 
00925     for (iloop = 0; iloop < MAXLOOP; iloop++) {
00926         n = 0;
00927 
00928         /* Ok, gather up all the stats */
00929 
00930         for (i = 0; i < nrows; i++) {
00931             
00932             /* Clipping criteria */
00933 
00934             if (ellipticity[i] < elllim && core_flux[i] < blim + 1.0 && 
00935                 core_flux[i] > flim + 1.0 && 
00936                 fabs(dc[i] - *medval) < 3.0*(*sigma) && cls[i] == -1.0)
00937                 work[n++] = dc[i];
00938         }
00939 
00940         /* Sort the work array and find the median and sigma */
00941 
00942         if (n > 2) {
00943             sort1(work,n);
00944             if (iloop == 0 && n > 10) {
00945                 anhist(work,n,medval,sigma);
00946             } else {
00947                 medstat(work,n,medval,sigma);
00948             }
00949         } else {
00950             *medval = 0.0;
00951             *sigma = 0.01;
00952         }
00953 
00954         /* Just in case... */
00955 
00956         *sigma = max(*sigma,0.01);
00957     }
00958 
00959     /* Tidy and exit */
00960 
00961     freespace(work);
00962     freespace(dc);
00963 }
00964 
00965 static void classstats_final(void) {
00966     int n,i,iloop,iarg,ii,iend,ncls,kk,k;
00967     float *work,ell,core,sig1,sig2,sig3,denom,w1,w2,w3,core_small;
00968     float core_large,*statistic,core_midd,pkht,xcor,cfit,csig;
00969     float *work1,junk,corlim1,corval1,corlim2,corval2,sigmaold;
00970 
00971     /* Initialise */
00972 
00973     sigma_final = 1.0e6;
00974     fit_final = 0.0;
00975     ncls = 0;
00976 
00977     /* Get some workspace */
00978 
00979     work = cpl_malloc(nrows*sizeof(float));
00980     work1 = cpl_malloc(nrows*sizeof(float));
00981     statistic = cpl_malloc(nrows*sizeof(float));
00982 
00983     /* Calculate the statistic now */
00984 
00985     for (i = 0; i < nrows; i++) {
00986         ell = ellipticity[i];
00987         pkht = peak_mag[i];
00988         core = core_flux[i];
00989         iarg = vircam_nint(10.0*(core - 5.0));
00990         iarg = max(1,min(NSAMPLE,iarg)) - 1;
00991         if (! poor) {
00992             sig1 = max(0.01,(fit1 - lower1[iarg])/3.0);
00993             sig2 = max(0.01,(fit2 - lower2[iarg])/3.0);
00994         } else {
00995             sig1 = max(0.01,(fit4 - lower1[iarg])/3.0);
00996             sig2 = max(0.01,(fit5 - lower2[iarg])/3.0);
00997         }
00998         sig3 = max(0.01,(fit3 - lower3[iarg])/3.0);
00999         denom = (wt1/sig1 + wt2/sig2 + wt3/sig3);
01000         w1 = (wt1/sig1)/denom;
01001         w2 = (wt2/sig2)/denom;
01002         w3 = (wt3/sig3)/denom;
01003         if (! poor) {
01004             core_small = core1_flux[i];
01005             core_large = core3_flux[i];
01006             statistic[i] = (core - core_small - fit1)*w1 + 
01007                 (core_large - core - fit2)*w2 + (core - pkht - fit3)*w3;
01008         } else {
01009             core_midd = core2_flux[i];
01010             core_large = core4_flux[i];
01011             statistic[i] = (core_midd - core - fit4)*w1 +
01012                 (core_large - core - fit5)*w2 + (core - pkht - fit3)*w3;
01013         }
01014     }
01015 
01016     /* Iteration loop.  Use only lower ellipticity images and relevant
01017        peak height range */
01018 
01019     for (iloop = 0; iloop < MAXLOOP; iloop++) {
01020         sigmaold = sigma_final;
01021         n = 0;
01022         for (i = 0; i < nrows ; i++) {
01023 
01024             ell = ellipticity[i];
01025             core = core_flux[i];
01026             if (ell < elllim && core < blim && core > flim && 
01027                 fabs((double)(statistic[i] - fit_final)) < 3.0*sigma_final) 
01028                 work[n++] = statistic[i];
01029 
01030             /* This information is to be used later to find the curvature of 
01031                saturated region */
01032 
01033             if (core > corlim && iloop == MAXLOOP-2) {
01034                 cls[ncls] = statistic[i];
01035                 sig[ncls++] = core;
01036             }
01037         }
01038 
01039         /* Median defines general fit */
01040 
01041         if (n > 2) {
01042             sort1(work,n);
01043             if (iloop == 0 && n > 10) {
01044                 anhist(work,n,&fit_final,&sigma_final);
01045             } else {
01046                 medstat(work,n,&fit_final,&sigma_final);
01047             }
01048             sigma_final = max(0.01,min(sigmaold,sigma_final));
01049         } else {
01050             fit_final = 0.0;
01051             sigma_final = 0.01;
01052         }
01053     }
01054 
01055     /* Now work out the curvature in the saturated region */
01056 
01057     sort2(sig,cls,ncls);
01058     ii = 0;
01059     xcor = 12.5;
01060     iend = 0;
01061     i = -1;
01062     corlim1 = 0.0;
01063     corlim2 = 0.0;
01064     corval1 = 0.0;
01065     corval2 = 0.0;
01066     while (iend == 0 && i < ncls-1) {
01067         i++;
01068         if (sig[i] > xcor+0.25 && ii >= 3) {
01069             medstat(work,ii,&cfit,&csig);
01070             for (iloop = 0; iloop < 3; iloop++) {
01071                 kk = 0;
01072                 for (k = 0; k < ii; k++) {
01073                     if (work[k] <= cfit + 3.0*csig)
01074                         work1[kk++] = work[k];
01075                 }
01076                 medstat(work1,kk,&cfit,&junk);
01077             }
01078             if (cfit <= fit_final + 3.0*sigma_final) {
01079                 corlim1 = xcor;
01080                 corval1 = cfit;
01081             } else {
01082                 corlim2 = xcor;
01083                 corval2 = cfit;
01084                 iend = 1;
01085             }
01086         } else {
01087             work[ii++] = cls[i];
01088         }
01089     }
01090 
01091     /* Estimate where core measure and statistic become unreliable */
01092 
01093     if (iend == 1) 
01094         corlim = corlim2 - 0.5*(corval2 - fit_final - 3.0*sigma_final)/(corval2 - corval1);
01095     else 
01096         corlim = corlim1;
01097     corlim = max(cormin,corlim);
01098     kk = 0;
01099     for (i = 0; i < nrows; i++) {
01100         core = core_flux[i];
01101         if (core >= corlim)
01102             work[kk++] = peak_height[i] + skylevel;
01103     }
01104     medstat(work,kk,&avsat,&junk);
01105     avsat = max(10000.0,avsat);
01106 
01107     /* Tidy and exit */
01108   
01109     freespace(work);
01110     freespace(work1);
01111     freespace(statistic);
01112 }
01113 
01114 
01115 static void medstat(float *array, int n, float *medval, float *sigval) {
01116     int lev1,lev2,lev3;
01117 
01118     /* Sort the array first, then choose the median.  The sigma is defined
01119        as half the distance between the two quartile points multiplied by
01120        the appropriate scaling factor (1.48) */
01121 
01122     if (n == 0) {
01123         *medval = 0.0;
01124         *sigval = 0.0;
01125         return;
01126     }
01127     sort1(array,n);
01128     lev1 = (int)(0.5*(float)(n + 1));
01129     lev2 = (int)(0.25*(float)(3*n + 3));
01130     lev3 = (int)(0.25*(float)(n + 3));
01131     *medval = array[lev1-1];
01132     *sigval = 1.48*0.5*(array[lev2-1] - array[lev3-1]);
01133 }
01134 
01135 
01136 static void sort1(float *a, int n) {
01137     int iii,ii,i,ifin,j;
01138     float b;
01139 
01140     iii = 2;
01141     while (iii < n)
01142         iii *= 2;
01143     iii = min(n,(3*iii)/4 - 1);
01144 
01145     while (iii > 1) {
01146         iii /= 2;
01147         ifin = n - iii;
01148         for (ii = 0; ii < ifin; ii++) {
01149             i = ii;
01150             j = i + iii;
01151             if (a[i] > a[j]) {
01152                 b = a[j];
01153                 while (1) {
01154                     a[j] = a[i];
01155                     j = i;
01156                     i = i - iii;
01157                     if (i < 0 || a[i] <= b) 
01158                         break;
01159                 }
01160                 a[j] = b;
01161             }
01162         }
01163     }
01164 }
01165 
01166 static void sort2(float *a1, float *a2, int n) {
01167     int iii,ii,i,ifin,j;
01168     float b1,b2;
01169 
01170     iii = 4;
01171     while (iii < n)
01172         iii *= 2;
01173     iii = min(n,(3*iii)/4 - 1);
01174 
01175     while (iii > 1) {
01176         iii /= 2;
01177         ifin = n - iii;
01178         for (ii = 0; ii < ifin; ii++) {
01179             i = ii;
01180             j = i + iii;
01181             if (a1[i] > a1[j]) {
01182                 b1 = a1[j];
01183                 b2 = a2[j];
01184                 while (1) {
01185                     a1[j] = a1[i];
01186                     a2[j] = a2[i];
01187                     j = i;
01188                     i = i - iii;
01189                     if (i < 0 || a1[i] <= b1) 
01190                         break;
01191                 }
01192                 a1[j] = b1;
01193                 a2[j] = b2;
01194             }
01195         }
01196     }
01197 }
01198 
01199 
01200 /*
01201 
01202 $Log: classify.c,v $
01203 Revision 1.9  2007/05/02 09:11:34  jim
01204 Modified to allow for inclusion of table WCS keywords into FITS header
01205 
01206 Revision 1.8  2007/03/01 12:38:26  jim
01207 Small modifications after a bit of code checking
01208 
01209 Revision 1.7  2006/06/13 14:06:57  jim
01210 The classification and statistic data rows must come from the original table
01211 rather than the copy as they are being written to
01212 
01213 Revision 1.6  2006/05/26 15:00:36  jim
01214 Now makes a copy of the input table so that it doesn't muck up the column
01215 values
01216 
01217 Revision 1.5  2006/05/18 12:34:20  jim
01218 Fixed header keywords for input information
01219 
01220 Revision 1.4  2006/03/15 10:43:42  jim
01221 Fixed a few things
01222 
01223 Revision 1.3  2005/11/28 13:50:06  jim
01224 Touched up a few things to make splint happy
01225 
01226 Revision 1.2  2005/11/03 13:28:51  jim
01227 All sorts of changes to tighten up error handling
01228 
01229 Revision 1.1  2005/09/22 08:41:03  jim
01230 first entry
01231 
01232 
01233 */

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