IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5084


Ignore:
Timestamp:
Sep 20, 2005, 6:48:16 PM (21 years ago)
Author:
eugene
Message:

alternative clipping for pmPSFtryMetric

Location:
trunk/psphot
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r5058 r5084  
    1515# LIBS  =       -lpsmodule $(LPSLIB)
    1616LIBS    =       $(LPSLIB)
    17 CFLAGS  =       $(INCS) -std=c99 -Wall -Werror -g
     17# CFLAGS        =       $(INCS) -std=c99 -Wall -Werror -g
     18CFLAGS  =       $(INCS) -std=c99 -g
    1819LFLAGS  =       $(LIBS)
    1920

  • trunk/psphot/src/pmPSFtry.c

    r5068 r5084  
    161161
    162162    // XXX this function wants aperture radius for pmSourcePhotometry
    163     pmPSFtryMetric (try, RADIUS);
     163    pmPSFtryMetric_Alt (try, RADIUS);
    164164    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    165165              modelName, try->psf->ApResid, try->psf->dApResid, try->psf->skyBias);
    … …  
    167167    return (try);
    168168}
    169 
    170169
    171170bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
    … …  
    195194  }
    196195  fclose (f);
     196
     197  // XXX EAM : try 3hi/1lo sigma clipping on the rflux v dap fit
    197198
    198199  // find min and max of (1/flux):
    … …  
    297298  return true;
    298299}
     300
     301bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS) {
     302
     303  // the measured (aperture - fit) magnitudes (dA == try->metric)
     304  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
     305  //     dA = dAo + dsky/flux
     306  //   where flux is the flux of the star
     307  // we fit this trend to find the infinite flux aperture correction (dAo),
     308  //   the nominal sky bias (dsky), and the error on dAo
     309  // the values of dA are contaminated by stars with close neighbors in the aperture
     310  //   we use an outlier rejection to avoid this bias
     311
     312  // rflux = ten(0.4*fitMag);
     313  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     314  for (int i = 0; i < try->sources->n; i++) {
     315    if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
     316    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
     317  }
     318
     319  // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
     320  psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     321
     322  // XXX EAM
     323  stats->min = 1.0;
     324  stats->max = 3.0;
     325  stats->clipIter = 3;
     326
     327  // linear fit to rfBin, daBin
     328  psPolynomial1D *poly = Polynomial1DAlloc_EAM (PS_POLYNOMIAL_ORD, 1);
     329
     330  // XXX EAM : this is the intended API (cycle 7? cycle 8?)
     331  poly = VectorClipFitPolynomial1D_EAM (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux);
     332
     333  fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     334
     335  // these stats should come back from the fit above...
     336  // psVector *daFit   = Polynomial1DEvalVector_EAM (poly, rflux);
     337  // psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) try->metric, "-", (void *) daFit);
     338
     339  // stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
     340  // stats->clipIter = 3;
     341  // stats->clipSigma = 3;
     342
     343  // stats = psVectorStats (stats, daResid, NULL, maskB, 1);
     344
     345  try->psf->ApResid = poly->coeff[0];
     346  try->psf->dApResid = stats->sampleStdev;
     347  try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     348
     349  FILE *f;
     350  f = fopen ("apresid.dat", "w");
     351  if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
     352
     353  for (int i = 0; i < try->sources->n; i++) {
     354    fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]);
     355  }
     356  fclose (f);
     357
     358  psFree (rflux);
     359  psFree (poly);
     360  psFree (stats);
     361
     362  // psFree (daFit);
     363  // psFree (daResid);
     364
     365  return true;
     366}

  • trunk/psphot/src/pmPSFtry.h

    r5058 r5084  
    2727pmPSFtry    *pmPSFtryModel (psArray *sources, char *modelName, float radius);
    2828bool         pmPSFtryMetric (pmPSFtry *try, float RADIUS);
     29bool         pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS);
    2930
    3031# endif

  • trunk/psphot/src/psPolynomials.c

    r5068 r5084  
    507507    psVector *zFit   = NULL;
    508508    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
    509 
    510     // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
     509    float minClipSigma = stats->min;
     510    float maxClipSigma = stats->max;
     511    float minClipValue;
     512    float maxClipValue;
     513
     514    // XXX EAM : use SAMPLE_MEDIAN and SAMPLE_STDEV for stats:
    511515    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    512516
    513     for (int N = 0; N < 3; N++) {
     517    for (int N = 0; N < stats->clipIter; N++) {
    514518        int Nkeep = 0;
    515519
    … …  
    519523
    520524        stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
    521         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
     525        minClipValue = -1*fabs(minClipSigma)*stats->sampleStdev;
     526        maxClipValue = +1*fabs(maxClipSigma)*stats->sampleStdev;
     527
     528        psTrace (".psphot.VectorClipFit", 4, "residual stats for robust fit:  %g +/- %g\n",
    522529                 stats->sampleMedian, stats->sampleStdev);
    523530
    … …  
    527534        for (int i = 0; i < zResid->n; i++) {
    528535            if (mask->data.U8[i]) continue;
    529             if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 2*stats->sampleStdev) {
     536            if (zResid->data.F64[i] - stats->sampleMedian > maxClipValue) {
     537                mask->data.U8[i] |= 0x01;
     538                continue;
     539            }       
     540            if (zResid->data.F64[i] - stats->sampleMedian < minClipValue) {
    530541                mask->data.U8[i] |= 0x01;
    531542                continue;
    … …  
    534545        }
    535546
    536         psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n",
     547        psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n",
    537548                 Nkeep, x->n);
    538549

  • trunk/psphot/src/psphot.h

    r5068 r5084  
    3939// output functions
    4040bool         pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    41 bool         pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    42 bool         pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    43 bool         pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    44 bool         pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    45 int          pmSourcesDophotType (pmSource *source);
     41bool         pmSourcesWriteOBJ  (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     42bool         pmSourcesWriteCMP  (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     43bool         pmSourcesWriteCMF  (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     44bool         pmSourcesWriteSX   (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     45
    4646bool         pmPeaksWriteText (psArray *sources, char *filename);
    4747bool         pmMomentsWriteText (psArray *sources, char *filename);
    … …  
    4949bool         pmModelWriteFLTs (psArray *sources, char *filename);
    5050bool         pmModelWriteNULLs (psArray *sources, char *filename);
     51
     52int          pmSourcesDophotType (pmSource *source);

  • trunk/psphot/src/psphotArguments.c

    r5049 r5084  
    3838  }
    3939
     40  // optional output residual image - add to config
     41  char *photcode = NULL;
     42  if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
     43    psArgumentRemove (N, argc, argv);
     44    photcode = psStringCopy (argv[N]);
     45    psArgumentRemove (N, argc, argv);
     46  }
     47
    4048  if (*argc != 4) usage ();
    4149
    … …  
    6068    psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
    6169  }
     70  if (photcode != NULL) {
     71    psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode);
     72  }
    6273  return (config);
    6374}
    … …  
    7081    fprintf (stderr, "  -weight (filename)\n");
    7182    fprintf (stderr, "  -resid (filename)\n");
     83    fprintf (stderr, "  -photcode (photcode)\n");
    7284    exit (2);
    7385}

  • trunk/psphot/src/psphotOutput.c

    r5073 r5084  
    2626
    2727    if (!strcasecmp (outputMode, "OBJ")) {
    28         pmSourcesWriteOBJ (imdata, outputFile, sources, psf, sky);
     28        pmSourcesWriteOBJ (imdata, config, outputFile, sources, psf, sky);
    2929        return;
    3030    }
    … …  
    3636 
    3737    if (!strcasecmp (outputMode, "CMP")) {
    38         pmSourcesWriteCMP (imdata, outputFile, sources, psf, sky);
     38        pmSourcesWriteCMP (imdata, config, outputFile, sources, psf, sky);
    3939        return;
    4040    }
    4141 
    4242    if (!strcasecmp (outputMode, "CMF")) {
    43         pmSourcesWriteCMF (imdata, outputFile, sources, psf, sky);
     43        pmSourcesWriteCMF (imdata, config, outputFile, sources, psf, sky);
    4444        return;
    4545    }
    … …  
    6868
    6969// dophot-style output list with fixed line width
    70 bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
     70bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    7171
    7272    int i, type, status;
    … …  
    7474    psF32 *PAR, *dPAR;
    7575    float sky, dmag, apMag, fitMag;
     76    float x, y, rflux;
     77    bool result;
    7678
    7779    psLine *line = psLineAlloc (104);  // 104 is dophot-defined line length
    … …  
    8284        return false;
    8385    }
     86
     87    float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS");
    8488
    8589    // write sources with models first
    … …  
    296300
    297301// elixir-style pseudo FITS table (header + ascii list)
    298 bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
     302bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    299303
    300304    int i, type, status;
    301305    pmModel *model;
     306    psMetadataItem *mdi;
    302307    psF32 *PAR, *dPAR;
    303     float sky, dmag, apMag, fitMag, lsky, Theta;
     308    float sky, dmag, apMag, fitMag, lsky;
     309    float x, y, rflux;
     310    bool result;
    304311
    305312    // create file, write-out header
    306313    unlink (filename);
    307314    psFits *fits = psFitsAlloc (filename);
     315
     316    // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0
     317    mdi = psMetadataLookup (imdata->header, "NAXIS");
     318    mdi->data.S32 = 0;
     319    mdi->type = PS_META_S32;
     320
     321    // psMetadataAdd (imdata->header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0);
    308322    psFitsWriteHeader (imdata->header, fits);
    309323    psFree (fits);
    … …  
    317331    fseek (f, 0, SEEK_END);
    318332
    319     psLine *line = psLineAlloc (66);  // 66 is imclean-defined line length
     333    psLine *line = psLineAlloc (67);  // 66 is imclean-defined line length
     334
     335    float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS");
     336
     337    // XXX EAM : need to add these concepts to the header (from psf and config)
     338    // fits_modify (header, "ZERO_PT", "%lf", 1, ZERO_POINT);
     339    // fits_modify (header, "FWHM_X", "%lf", 1, FWHMx);
     340    // fits_modify (header, "FWHM_Y", "%lf", 1, FWHMy);
     341    // fits_modify (header, "ANGLE", "%lf", 1, angle);
     342    // fits_modify (header, "FSATUR", "%lf", 1, saturate);  // can use 0.0
     343    // fits_modify (header, "FLIMIT", "%lf", 1, complete);  // can use 0.0
     344    // fits_modify (header, "NSTARS", "%d", 1, N);
    320345
    321346    // write sources with models first
    … …  
    404429        psLineAdd (line, "%6.2f ", PAR[4]);  // should be 'FHWM x'
    405430        psLineAdd (line, "%6.2f ", PAR[5]);  // should be 'FHWM y'
    406         psLineAdd (line, "%5.1f\n", Theta);   // should be theta
     431        psLineAdd (line, "%5.1f\n", 0);   // should be theta
    407432        fwrite (line->line, 1, line->Nline, f);
    408433    }
    … …  
    414439// this format consists of a header derived from the image header
    415440// followed by a zero-size matrix, followed by the table data
    416 bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
     441bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    417442
    418443    psArray *table;
     444    psMetadataItem *mdi;
    419445    psMetadata *row;
    420446    psMetadata *theader;
     447    int i, type, status;
     448    pmModel *model;
     449    psF32 *PAR, *dPAR;
     450    float sky, dmag, apMag, fitMag;
     451    float x, y, rflux, lsky;
     452    bool result;
    421453
    422454    // set NAXIS to 0
    423     psMetadataAdd (header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0);
     455    fprintf (stderr, "setting naxis\n");
     456    mdi = psMetadataLookup (imdata->header, "NAXIS");
     457    mdi->data.S32 = 0;
     458    mdi->type = PS_META_S32;
     459
     460    // psMetadataAdd (imdata->header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0);
    424461
    425462    table = psArrayAlloc (sources->n);
    426463    sources->n = 0;
     464
     465    float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS");
    427466
    428467    for (i = 0; i < sources->n; i++) {
    … …  
    499538        lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]);
    500539
    501         table = psMetdataAlloc ();
    502         psMetadataAdd (table, PS_LIST_HEAD, "X_PIX",   PS_META_F32, "", PAR[2]);
    503         psMetadataAdd (table, PS_LIST_HEAD, "Y_PIX",   PS_META_F32, "", PAR[3]);
    504         psMetadataAdd (table, PS_LIST_HEAD, "MAG_RAW", PS_META_F32, "", fitMag + 25.0);
    505         psMetadataAdd (table, PS_LIST_HEAD, "MAG_ERR", PS_META_F32, "", (int)(1000*dmag));
    506         psMetadataAdd (table, PS_LIST_HEAD, "MAG_GAL", PS_META_F32, "", type);
    507         psMetadataAdd (table, PS_LIST_HEAD, "MAG_AP",  PS_META_F32, "", lsky);
    508         psMetadataAdd (table, PS_LIST_HEAD, "LOG_SKY", PS_META_F32, "", 32.0);    // should be 'Mgal'
    509         psMetadataAdd (table, PS_LIST_HEAD, "FWHM_X",  PS_META_F32, "", apMag + 25.0); // 25.0 should come from config
    510         psMetadataAdd (table, PS_LIST_HEAD, "FWHM_Y",  PS_META_F32, "", PAR[4]);  // should be 'FHWM x'
    511         psMetadataAdd (table, PS_LIST_HEAD, "THETA",   PS_META_F32, "", PAR[5]);  // should be 'FHWM y'
    512         psMetadataAdd (table, PS_LIST_HEAD, "DOPHOT",  PS_META_S8,  "", Theta);   // should be theta
     540        row = psMetadataAlloc ();
     541        psMetadataAdd (row, PS_LIST_HEAD, "X_PIX",   PS_META_F32, "", PAR[2]);
     542        psMetadataAdd (row, PS_LIST_HEAD, "Y_PIX",   PS_META_F32, "", PAR[3]);
     543        psMetadataAdd (row, PS_LIST_HEAD, "MAG_RAW", PS_META_F32, "", fitMag + 25.0);
     544        psMetadataAdd (row, PS_LIST_HEAD, "MAG_ERR", PS_META_F32, "", (int)(1000*dmag));
     545        psMetadataAdd (row, PS_LIST_HEAD, "MAG_GAL", PS_META_F32, "", type);
     546        psMetadataAdd (row, PS_LIST_HEAD, "MAG_AP",  PS_META_F32, "", lsky);
     547        psMetadataAdd (row, PS_LIST_HEAD, "LOG_SKY", PS_META_F32, "", 32.0);    // should be 'Mgal'
     548        psMetadataAdd (row, PS_LIST_HEAD, "FWHM_X",  PS_META_F32, "", apMag + 25.0); // 25.0 should come from config
     549        psMetadataAdd (row, PS_LIST_HEAD, "FWHM_Y",  PS_META_F32, "", PAR[4]);  // should be 'FHWM x'
     550        psMetadataAdd (row, PS_LIST_HEAD, "THETA",   PS_META_F32, "", PAR[5]);  // should be 'FHWM y'
     551        psMetadataAdd (row, PS_LIST_HEAD, "DOPHOT",  PS_TYPE_S8,  "", 0);   // should be theta
    513552        // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY",   PS_META_STR, "", NULL);
    514553   
    … …  
    517556    }
    518557
    519     image = psImageAlloc (0, 0, PS_TYPE_F32);
    520     theader = psFitsHeaderFromTable (NULL, table, "SMPFILE");
     558    // psAbort ("psphotOutput", "FITS table output not finished");
     559
     560    psImage *image = psImageAlloc (1, 1, PS_TYPE_F32);
     561    // XXX missing: theader = psFitsHeaderFromTable (NULL, table, "SMPFILE");
     562
     563    theader = psMetadataAlloc ();
     564    psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_META_STR, "extension name", "SMPFILE");
    521565
    522566    unlink (filename);
    523567    psFits *fits = psFitsAlloc (filename);
    524     psFitsWriteImage (fits, header, image, 0, NULL);
     568    psFitsWriteImage (fits, imdata->header, image, 0);
    525569    psFitsWriteTable (fits, theader, table);
    526570    psFree (fits);
Note: See TracChangeset for help on using the changeset viewer.