00001
00002
00003
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
00031
00032 static long nrows;
00033 static float thresh,skylevel,skynoise,rcore,exptime;
00034
00035
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
00046
00047 static int nstar,ngal,njunk,ncmp;
00048
00049
00050
00051 static float avsat,corlim,cormin,apcpkht,apcor,apcor1,apcor2,apcor3,apcor4;
00052 static float apcor5,apcor6,apcor7;
00053
00054
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
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
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
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
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
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
00142
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
00160
00161 workspace = cpl_malloc(2*nrows*sizeof(float));
00162 peak_mag = workspace;
00163 work = workspace + nrows;
00164
00165
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
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
00187
00188 poor = 0;
00189 if (fwhm > max(5.0,rcore*sqrt(2.0)))
00190 poor = 1;
00191
00192
00193
00194 classify_run();
00195
00196
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
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
00228
00229 freespace(workspace);
00230 freetable(catcopy);
00231
00232
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
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
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
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
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
00326
00327 histo = cpl_calloc(MAXHIST,sizeof(int));
00328 sval = cpl_calloc(MAXHIST,sizeof(float));
00329
00330
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
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
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
00365
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
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
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
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
00420
00421 work = cpl_malloc(nrows*sizeof(float));
00422
00423
00424
00425 lower[0] = cmin;
00426 lower[1] = cmax;
00427 asign = ((small == 1) ? -1.0 : 1.0);
00428
00429
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
00448
00449 medstat(work,n,avsig,&junk);
00450 freespace(work);
00451
00452
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
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
00486
00487 work = cpl_malloc(nrows*sizeof(float));
00488
00489
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
00500
00501 medstat(work,n,avsig,&junk);
00502 freespace(work);
00503 *wt = min(5.0,max(1.0,*avsig/sigma));
00504
00505
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
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
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
00552
00553 classstats_el();
00554
00555
00556
00557
00558 classstats(core_flux,core1_flux,1,&fit1,&sigma1);
00559
00560
00561
00562 classstats(core_flux,core3_flux,0,&fit2,&sigma2);
00563
00564
00565
00566 classstats(core_flux,core2_flux,0,&fit4,&sigma4);
00567
00568
00569
00570 classstats(core_flux,core4_flux,0,&fit5,&sigma5);
00571
00572
00573
00574 classstats(core_flux,peak_mag,1,&fit3,&sigma3);
00575
00576
00577
00578 classstats_ellf(fluxlim);
00579
00580
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
00590
00591
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
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
00604
00605 boundpk(core_flux,peak_mag,fit3,sigma3,&wt3,&avsig3,lower3,upper3);
00606
00607
00608
00609
00610 classstats_final();
00611
00612
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
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
00680
00681 sigtot = (fit_final - lower[iarg])/5.0;
00682 sig[i] = (statistic - fit_final)/sigtot;
00683
00684
00685
00686
00687 if (core - pkht - fit3 < -4.0*sig3)
00688 cls[i] = 0.0;
00689
00690
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
00697
00698 if (core > corlim && statistic >= lower[iarg])
00699 cls[i] = -1.0;
00700
00701
00702
00703 if (ell > 0.9 && core < corlim)
00704 cls[i] = 0.0;
00705
00706
00707
00708 if (core < fluxlim)
00709 cls[i] = 0.0;
00710
00711
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
00724
00725 classstats_ap(&fit0,&sigma0);
00726 fit0 = max(fit5,fit0);
00727 apcpkht = fit0 + fit3;
00728 switch (ncols) {
00729 case 32:
00730 apcor1 = fit0 + fit1;
00731 apcor = fit0;
00732 apcor2 = fit0 - fit4;
00733 apcor3 = fit0 - fit2;
00734 apcor4 = fit0 - fit5;
00735 apcor5 = 0.0;
00736 break;
00737 case 80:
00738 apcor1 = fit0 + fit1;
00739 apcor2 = fit0 + 0.5*fit1;
00740 apcor3 = fit0;
00741 apcor4 = fit0 - fit4;
00742 apcor5 = fit0 - fit2;
00743 apcor6 = fit0 - fit5;
00744 apcor7 = 0.0;
00745 break;
00746 }
00747
00748
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
00769
00770 *medval = 0.0;
00771 *sigma = 1.0e6;
00772 amult = (small == 1 ? -1.0 : 1.0);
00773
00774
00775
00776 work = cpl_malloc(nrows*sizeof(float));
00777 dc = cpl_malloc(nrows*sizeof(float));
00778
00779
00780
00781 for (i = 0; i < nrows; i++)
00782 dc[i] = amult*(core2[i] - core1[i]);
00783
00784
00785
00786 for (iloop = 0; iloop < MAXLOOP; iloop++) {
00787 sigmaold = *sigma;
00788 n = 0;
00789
00790
00791
00792 for (i = 0; i < nrows; i++) {
00793
00794
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
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
00818
00819 *sigma = max(*sigma,0.01);
00820 }
00821
00822
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
00833
00834 sigell = 1.0e6;
00835 fitell = 0.0;
00836
00837
00838
00839 work = cpl_malloc(nrows*sizeof(float));
00840
00841
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
00860
00861 freespace(work);
00862 }
00863
00864 static void classstats_ellf(float fluxlim) {
00865 int iloop,n,i;
00866 float *work;
00867
00868
00869
00870 sigellf = 1.0e6;
00871 fitellf = 0.0;
00872
00873
00874
00875 work = cpl_malloc(nrows*sizeof(float));
00876
00877
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
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
00906
00907 *medval = 0.0;
00908 *sigma = 1.0e6;
00909 elllim = min(0.5,max(0.2,fitell+2.0*sigell));
00910
00911
00912
00913 work = cpl_malloc(nrows*sizeof(float));
00914 dc = cpl_malloc(nrows*sizeof(float));
00915
00916
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
00924
00925 for (iloop = 0; iloop < MAXLOOP; iloop++) {
00926 n = 0;
00927
00928
00929
00930 for (i = 0; i < nrows; i++) {
00931
00932
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
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
00955
00956 *sigma = max(*sigma,0.01);
00957 }
00958
00959
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
00972
00973 sigma_final = 1.0e6;
00974 fit_final = 0.0;
00975 ncls = 0;
00976
00977
00978
00979 work = cpl_malloc(nrows*sizeof(float));
00980 work1 = cpl_malloc(nrows*sizeof(float));
00981 statistic = cpl_malloc(nrows*sizeof(float));
00982
00983
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
01017
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
01031
01032
01033 if (core > corlim && iloop == MAXLOOP-2) {
01034 cls[ncls] = statistic[i];
01035 sig[ncls++] = core;
01036 }
01037 }
01038
01039
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
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
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
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
01119
01120
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
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233