IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5068


Ignore:
Timestamp:
Sep 19, 2005, 12:26:40 PM (21 years ago)
Author:
eugene
Message:

fixing errors related to apResid and skyBias

Location:
trunk/psphot/src
Files:
1 deleted
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/pmPSF.c

    r5058 r5068  
    4040        // XXX EAM : make this a user-defined value?
    4141        // XXX EAM : future version (0.7.0?) psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
    42         psf->params->data[i] = Polynomial2DAlloc_EAM(PS_POLYNOMIAL_ORD, 1, 1);
     42        psf->params->data[i] = Polynomial2DAlloc_EAM(PS_POLYNOMIAL_ORD, 0, 0);
    4343    }
    4444

  • trunk/psphot/src/pmPSFtry.c

    r5058 r5068  
    183183  //   we use an outlier rejection to avoid this bias
    184184
     185  FILE *f;
     186  f = fopen ("apresid.dat", "w");
     187  if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
     188
    185189  // rflux = ten(0.4*fitMag);
    186190  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
    … …  
    188192    if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    189193    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
     194    fprintf (f, "%3d %8.4f %12.5e %8.4f\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i]);
    190195  }
     196  fclose (f);
    191197
    192198  // find min and max of (1/flux):
    … …  
    236242    // measure statistics only on upper 50% of points
    237243    // this would be easier if we could sort in reverse:
    238     //
    239     // psVectorSort (tmp, tmp);
    240     // tmp->n = 0.5*tmp->n;
    241     // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    242     // psVectorStats (stats, tmp, NULL, NULL, 0);
    243     // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    244244
    245245    psVectorSort (tmp, tmp);
    … …  
    267267
    268268  // XXX EAM : this is the intended API (cycle 7? cycle 8?)
    269   poly = VectorFitPolynomial1D_EAM (poly, maskB, 1, daBin, NULL, rfBin);
     269  psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     270  poly = VectorClipFitPolynomial1D_EAM (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
     271
     272  // poly = VectorFitPolynomial1D_EAM (poly, maskB, 1, daBin, NULL, rfBin);
    270273
    271274  // XXX EAM : replace this when the above version is implemented
    … …  
    290293  psFree (poly);
    291294  psFree (stats);
     295  psFree (fitstat);
    292296
    293297  return true;

  • trunk/psphot/src/psLibUtils.h

    r5058 r5068  
    6262psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
    6363psPolynomial1D *VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly, psVector* mask, psMaskType maskValue, const psVector* y, const psVector* yErr, const psVector* x);
     64psPolynomial1D *VectorClipFitPolynomial1D_EAM(psPolynomial1D* poly, psStats *stats, psVector* mask, psMaskType maskValue, const psVector* z, const psVector* zErr, const psVector* x);
    6465
    6566psPolynomial2D *Polynomial2DAlloc_EAM(psPolynomialType type, psS32 nXorder, psS32 nYorder);

  • trunk/psphot/src/psPolynomials.c

    r4977 r5068  
    497497
    498498// XXX EAM : be careful here with F32 vs F64 vectors
     499psPolynomial1D* VectorClipFitPolynomial1D_EAM(psPolynomial1D* poly,
     500                                              psStats *stats,
     501                                              psVector* mask,
     502                                              psMaskType maskValue,
     503                                              const psVector* z,
     504                                              const psVector* dz,
     505                                              const psVector* x)
     506{
     507    psVector *zFit   = NULL;
     508    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
     509
     510    // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
     511    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     512
     513    for (int N = 0; N < 3; N++) {
     514        int Nkeep = 0;
     515
     516        poly   = VectorFitPolynomial1D_EAM (poly, mask, maskValue, z, dz, x);
     517        zFit   = Polynomial1DEvalVector_EAM (poly, x);
     518        zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
     519
     520        stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
     521        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
     522                 stats->sampleMedian, stats->sampleStdev);
     523
     524        // set mask if pts are not valid
     525        // we are masking out any point which is out of range
     526        // recovery is not allowed with this scheme
     527        for (int i = 0; i < zResid->n; i++) {
     528            if (mask->data.U8[i]) continue;
     529            if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 2*stats->sampleStdev) {
     530                mask->data.U8[i] |= 0x01;
     531                continue;
     532            }       
     533            Nkeep ++;
     534        }
     535
     536        psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n",
     537                 Nkeep, x->n);
     538
     539        psFree (zFit);
     540    }
     541    psFree (zResid);
     542    return (poly);
     543}
     544
     545
     546// XXX EAM : be careful here with F32 vs F64 vectors
    499547psPolynomial2D* VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly,
    500548                                              psStats *stats,
    … …  
    510558
    511559    // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
    512     stats->options |= (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     560    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    513561
    514562    for (int N = 0; N < 3; N++) {
     563        int Nkeep = 0;
     564
    515565        poly   = VectorFitPolynomial2D_EAM (poly, mask, maskValue, z, dz, x, y);
    516566        zFit   = Polynomial2DEvalVectorD_EAM (poly, x, y);
    … …  
    519569        stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
    520570        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
    521                  stats->sampleMean, stats->sampleStdev);
     571                 stats->sampleMedian, stats->sampleStdev);
    522572
    523573        // set mask if pts are not valid
    … …  
    526576        for (int i = 0; i < zResid->n; i++) {
    527577            if (mask->data.U8[i]) continue;
    528             if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
    529                 mask->data.U8[i] &= 0x01;
     578            if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 3*stats->sampleStdev) {
     579                mask->data.U8[i] |= 0x01;
    530580                continue;
    531581            }       
    532         }
     582            Nkeep ++;
     583        }
     584
     585        psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n",
     586                 Nkeep, x->n);
     587
     588//      psTrace (".psphot.RobustFit", 4, "model pars: %f %f %f %f\n",
     589//               poly->coeff[0][0], poly->coeff[1][0],
     590//               poly->coeff[0][1], poly->coeff[1][1]);
     591
     592        psTrace (".psphot.RobustFit", 4, "model pars: %f\n", poly->coeff[0][0]);
     593
    533594        psFree (zFit);
    534595    }

  • trunk/psphot/src/psphot.h

    r5058 r5068  
    4242bool         pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    4343bool         pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    44 bool         pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     44bool         pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    4545int          pmSourcesDophotType (pmSource *source);
    4646bool         pmPeaksWriteText (psArray *sources, char *filename);

  • trunk/psphot/src/psphotOutput.c

    r5058 r5068  
    3131 
    3232    if (!strcasecmp (outputMode, "SX")) {
    33         pmSourcesWriteSX (imdata, outputFile, sources, psf, sky);
     33        pmSourcesWriteSX (imdata, config, outputFile, sources, psf, sky);
    3434        return;
    3535    }
    … …  
    159159}
    160160
    161 // elixir/sextractor-style output list with fixed line width
    162 bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
     161// elixir-mode / sextractor-style output list with fixed line width
     162bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    163163
    164164    int i, status, flags;
    … …  
    166166    psF32 *PAR, *dPAR;
    167167    float sky, dmag, apMag, fitMag;
    168 
    169     psLine *line = psLineAlloc (110);  // 104 is dophot-defined line length
     168    float x, y, rflux;
     169    bool result;
     170
     171    psLine *line = psLineAlloc (110);  // 110 is sextractor line length
    170172
    171173    FILE *f = fopen (filename, "w");
    … …  
    174176        return false;
    175177    }
     178
     179    // aperture radius for ap magnitude
     180    float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS");
    176181
    177182    // write sources with models first
    … …  
    214219        dPAR = model->dparams->data.F32;
    215220
     221        x = source->peak->x;
     222        y = source->peak->y;
     223
     224        // we have already (psphotApplyPSF) defined pixels at least OUTER_RADIUS from source
     225        // we need to mask pixels to measure the aperture magnitude
     226        // set aperture mask circle of PSF_FIT_RADIUS
     227        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
     228
    216229        // save local sky for later
    217230        sky = model->params->data.F32[0];
    … …  
    224237
    225238        // measure object photometry
    226         status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
     239        status  = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
     240
     241        // correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS
     242        rflux   = pow (10.0, 0.4*fitMag);
     243        apMag  -= rflux * psf->skyBias * (M_PI * PS_SQR(RADIUS));
    227244        fitMag += psf->ApResid;
    228245
    … …  
    232249        pmSourceSubModel (source->pixels, source->mask, model, false);
    233250        model->params->data.F32[0] = sky;
     251
     252        // unmask aperture
     253        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
    234254
    235255        if (status == FALSE) continue;
Note: See TracChangeset for help on using the changeset viewer.