IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4375


Ignore:
Timestamp:
Jun 24, 2005, 5:38:37 AM (21 years ago)
Author:
eugene
Message:

adding model tests

Location:
trunk/psphot
Files:
3 added
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r4251 r4375  
    4141$(SRC)/find_defects.$(ARCH).o
    4242
     43FITSOURCE = \
     44$(SRC)/fitsource.$(ARCH).o \
     45$(SRC)/onesource.$(ARCH).o \
     46$(SRC)/fs_args.$(ARCH).o \
     47$(SRC)/psphot-utils.$(ARCH).o \
     48$(SRC)/psPolynomials.$(ARCH).o \
     49$(SRC)/psUtils.$(ARCH).o \
     50$(SRC)/setup.$(ARCH).o \
     51$(SRC)/LocalSky.$(ARCH).o
     52
    4353psphot: $(BIN)/psphot.$(ARCH)
     54$(BIN)/psphot.$(ARCH) : $(PSPHOT)
     55$(PSPHOT): $(SRC)/psphot.h
    4456
    45 $(BIN)/psphot.$(ARCH) : $(PSPHOT)
    46 
    47 $(PSPHOT): $(SRC)/psphot.h
     57fitsource: $(BIN)/fitsource.$(ARCH)
     58$(BIN)/fitsource.$(ARCH) : $(FITSOURCE)
     59$(FITSOURCE): $(SRC)/psphot.h
    4860
    4961INSTALL = psphot

  • trunk/psphot/src/LocalSky.c

    r4251 r4375  
    1313    srcRegion = psRegionForImage (srcRegion, imdata->image, srcRegion);
    1414
    15 //  psImage *subImage = psImageSubset(imdata->image, srcRegion);
    16 //  psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
    17 //  psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
     15    //  use these when psImageSubset is updated
     16    //  psImage *subImage = psImageSubset(imdata->image, srcRegion);
     17    //  psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
     18    //  psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
    1819
    1920    psImage *subImage = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
    … …  
    4748    psImageMaskRegion (mask, srcRegion, AND, 0x7f);
    4849
    49     source->moments = pmMomentsAlloc();
    5050    psF64 tmpF64;
    5151    p_psGetStatValue(myStats, &tmpF64);
     52    if (isnan(tmpF64)) return (false);
     53
     54    source->moments = pmMomentsAlloc();
    5255    source->moments->Sky = (psF32) tmpF64;
    5356

  • trunk/psphot/src/apply_psf_model.c

    r4251 r4375  
    1111    float y;
    1212    int   Nfit = 0;
     13    int   Nsub = 0;
    1314    int   Niter = 0;
    1415
    … …  
    7172
    7273        mark_psf_source (source, shapeNsigma);
    73         subtract_psf_source (source);
     74        if (subtract_psf_source (source)) {
     75          Nsub ++;
     76        }
    7477    }
    7578    psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
     79    psLogMsg ("psphot", 3, "subtracted %d PSF objects\n", Nsub);
    7680    return (true);
    7781}

  • trunk/psphot/src/by_SN.c

    r4129 r4375  
    99    psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
    1010    psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN;
     11    if (isnan (fA)) fA = 0;
     12    if (isnan (fB)) fB = 0;
    1113
    1214    psF32 diff = fA - fB;

  • trunk/psphot/src/choose_psf_model.c

    r4251 r4375  
    33// test PSF models and select best option
    44
    5 pmPSF *choose_psf_model (psMetadata *config, psArray *sources)
     5pmPSF *choose_psf_model (psMetadata *config, psArray *sources, psStats *skystats)
    66{
    77    bool        status;
    … …  
    3434        modelName = psMetadataLookupPtr (&status, config, key);
    3535        tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS);
     36
     37        // dump psf test data to file:
     38        {
     39          sprintf (key, "metric.%d.dat", i);
     40          FILE *f = fopen (key, "w");
     41          test = tests->data[i];
     42
     43
     44          for (int j = 0; j < test->sources->n; j++) {
     45            psSource *source = test->sources->data[j];
     46            psModel  *model  = test->modelPSF->data[j];     
     47            float sky = 0;
     48
     49            if (model != NULL) sky = model->params->data.F32[0];
     50
     51            fprintf (f, "%3d %6.1f %6.1f  %7.4f %7.4f  %8.4f\n",
     52                     j, source->moments->x, source->moments->y,
     53                     test->metric->data.F64[j], test->fitMag->data.F64[j], sky);
     54          }
     55          fclose (f);
     56        }
    3657    }
    3758
    … …  
    4061    test = tests->data[0];
    4162    int   bestN = 0;
    42     float bestM = test->metricStats->clippedStdev;
     63    float bestM = test->dApResid;
    4364    for (int i = 1; i < Ntest; i++) {
    4465        test = tests->data[i];
    45         int M = test->metricStats->clippedStdev;
     66        int M = test->dApResid;
    4667        if (M < bestM) {
    4768            bestM = M;
    … …  
    5172    test = tests->data[bestN];
    5273    modelName = psModelGetType (test->modelType);
    53     psLogMsg ("psphot.pspsf", 3, "selected psf model %s, metric: %f +/- %f\n", modelName, test->metricStats->clippedMean, test->metricStats->clippedStdev);
     74    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid);
    5475
    5576    // set source->model based on best psf model fit

  • trunk/psphot/src/load_args.c

    r4251 r4375  
    4545    fprintf (stderr, "USAGE: psphot (image.fits) (output.fits) (config)\n");
    4646    fprintf (stderr, "options: \n");
    47     fprintf (stderr, "  --mask  (filename)");
    48     fprintf (stderr, "  --noise (filename)");
     47    fprintf (stderr, "  --mask  (filename)\n");
     48    fprintf (stderr, "  --noise (filename)\n");
    4949    exit (2);
    5050}

  • trunk/psphot/src/mark_psf_source.c

    r4129 r4375  
    1414// any object which meets the criterion is marked as
    1515// PS_SOURCE_BRIGHTSTAR
     16# define MIN_DS 0.01
    1617bool mark_psf_source (psSource *source, float shapeNsigma)
    1718{
    … …  
    3031    dSY = source->modelPSF->dparams->data.F32[5];
    3132
    32     nSx = dSX * SX * SN;
    33     nSy = dSY * SY * SN;
     33    nSx = dSX / hypot (MIN_DS, 1 / (SX * SN));
     34    nSy = dSY / hypot (MIN_DS, 1 / (SY * SN));
     35    // sigma = 1 / SX
     36    // dsx = 1 / (SX * SN)
     37    // dsx_o = hypot (1/SX*SN + MIN_DSX)
    3438
    3539    // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]

  • trunk/psphot/src/psPolynomials.c

    r4129 r4375  
    11# include "psphot.h"
    22
    3 // XXX EAM : this version uses myPoly->nX as Norder, not Nterms
    4 psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
    5                                  const psVector *x,
    6                                  const psVector *y)
    7 
    8 {
    9     PS_POLY_CHECK_NULL(myPoly, NULL);
    10     PS_VECTOR_CHECK_NULL(x, NULL);
    11     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
    12     PS_VECTOR_CHECK_NULL(y, NULL);
    13     PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
    14 
    15     psVector *tmp;
    16     psS32 vecLen=x->n;
    17 
    18     // Determine the length of the output vector to by the minimum of the x,y vectors
    19     if (y->n < vecLen) {
    20         vecLen = y->n;
    21     }
    22 
    23     // Create output vector to return
    24     tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
    25 
    26     // Evaluate the polynomial at the specified points
    27     for (psS32 i=0; i<vecLen; i++) {
    28         tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
    29     }
    30 
    31     // Return output vector
    32     return(tmp);
    33 }
    34 
    35 // XXX EAM : this version uses the F64 vectors
    36 psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
    37                                   const psVector *x,
    38                                   const psVector *y)
    39 
     3// write out the terms of the given 1D polynomial
     4void psPolynomial1DDump (psPolynomial1D *poly) {
     5
     6    for (int i = 0; i < poly->n + 1; i++) {
     7        fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]);
     8    }
     9}   
     10
     11psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly)
     12{
     13    psS32 loop_x = 0;
     14    psF32 polySum = 0.0;
     15    psF32 xSum = 1.0;
     16
     17    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
     18        if (myPoly->mask[loop_x] == 0) {
     19            polySum += xSum * myPoly->coeff[loop_x];
     20        }
     21        xSum *= x;
     22    }
     23
     24    return(polySum);
     25}
     26
     27psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
     28                                       const psVector *x)
    4029{
    4130    PS_POLY_CHECK_NULL(myPoly, NULL);
    4231    PS_VECTOR_CHECK_NULL(x, NULL);
    4332    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
    44     PS_VECTOR_CHECK_NULL(y, NULL);
    45     PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
    4633
    4734    psVector *tmp;
    48     psS32 vecLen=x->n;
    49 
    50     // Determine the length of the output vector to by the minimum of the x,y vectors
    51     if (y->n < vecLen) {
    52         vecLen = y->n;
    53     }
    54 
    55     // Create output vector to return
    56     tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    57 
    58     // Evaluate the polynomial at the specified points
    59     for (psS32 i=0; i<vecLen; i++) {
    60         tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
    61     }
    62 
    63     // Return output vector
     35
     36    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
     37    for (psS32 i=0;i<x->n;i++) {
     38        tmp->data.F64[i] = Polynomial1DEval_EAM(x->data.F64[i], myPoly);
     39    }
     40
    6441    return(tmp);
    6542}
    … …  
    9370}
    9471
    95 // XXX EAM : use Nterm = Norder + 1 definition 
    96 // the user requests a polynomial of order Norder
    97 psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
    98                                     psPolynomialType type)
    99 {
    100     PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
    101     PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
    102 
    103     psS32 x = 0;
    104     psS32 y = 0;
    105     psS32 nXterm = nXorder + 1;
    106     psS32 nYterm = nYorder + 1;
    107     psPolynomial2D* newPoly = NULL;
    108 
    109     newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
    110     // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
    111     // XXX EAM : me, being lazy
    112 
    113     newPoly->type = type;
    114     newPoly->nX = nXorder;
    115     newPoly->nY = nYorder;
    116 
    117     newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
    118     newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
    119     newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
    120     for (x = 0; x < nXterm; x++) {
    121         newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
    122         newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
    123         newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
    124     }
    125     for (x = 0; x < nXterm; x++) {
    126         for (y = 0; y < nYterm; y++) {
    127             newPoly->coeff[x][y] = 0.0;
    128             newPoly->coeffErr[x][y] = 0.0;
    129             newPoly->mask[x][y] = 0;
    130         }
    131     }
    132     return(newPoly);
    133 }
    134 
    13572// XXX EAM : my alternate BuildSums1D
    13673static psVector *BuildSums1D(psVector* sums,
    … …  
    15390        sums->data.F64[i] = xSum;
    15491        xSum *= x;
    155     }
    156     return (sums);
    157 }
    158 
    159 // XXX EAM : BuildSums2D in analogy with BuildSums1D
    160 static psImage *BuildSums2D(psImage* sums,
    161                             psF64 x,      psF64 y,
    162                             psS32 nXterm, psS32 nYterm)
    163 {
    164     psS32 nXsum = 0;
    165     psS32 nYsum = 0;
    166     psF64 xSum = 1.0;
    167     psF64 ySum = 1.0;
    168 
    169     nXsum = 2*nXterm;
    170     nYsum = 2*nYterm;
    171     if (sums == NULL) {
    172         sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
    173     }
    174     if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
    175         psFree (sums);
    176         sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
    177     }
    178 
    179     ySum = 1.0;
    180     for (int j = 0; j < nYsum; j++) {
    181         xSum = ySum;
    182         for (int i = 0; i < nXsum; i++) {
    183             sums->data.F64[i][j] = xSum;
    184             xSum *= x;
    185         }
    186         ySum *= y;
    18792    }
    18893    return (sums);
    … …  
    291196}
    292197
     198// ********************** 2D polynomial functions ******************
     199
     200// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
     201psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
     202                                 const psVector *x,
     203                                 const psVector *y)
     204
     205{
     206    PS_POLY_CHECK_NULL(myPoly, NULL);
     207    PS_VECTOR_CHECK_NULL(x, NULL);
     208    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     209    PS_VECTOR_CHECK_NULL(y, NULL);
     210    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     211
     212    psVector *tmp;
     213    psS32 vecLen=x->n;
     214
     215    // Determine the length of the output vector to by the minimum of the x,y vectors
     216    if (y->n < vecLen) {
     217        vecLen = y->n;
     218    }
     219
     220    // Create output vector to return
     221    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
     222
     223    // Evaluate the polynomial at the specified points
     224    for (psS32 i=0; i<vecLen; i++) {
     225        tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
     226    }
     227
     228    // Return output vector
     229    return(tmp);
     230}
     231
     232// XXX EAM : this version uses the F64 vectors
     233psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
     234                                  const psVector *x,
     235                                  const psVector *y)
     236
     237{
     238    PS_POLY_CHECK_NULL(myPoly, NULL);
     239    PS_VECTOR_CHECK_NULL(x, NULL);
     240    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     241    PS_VECTOR_CHECK_NULL(y, NULL);
     242    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     243
     244    psVector *tmp;
     245    psS32 vecLen=x->n;
     246
     247    // Determine the length of the output vector to by the minimum of the x,y vectors
     248    if (y->n < vecLen) {
     249        vecLen = y->n;
     250    }
     251
     252    // Create output vector to return
     253    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
     254
     255    // Evaluate the polynomial at the specified points
     256    for (psS32 i=0; i<vecLen; i++) {
     257        tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
     258    }
     259
     260    // Return output vector
     261    return(tmp);
     262}
     263
     264// XXX EAM : use Nterm = Norder + 1 definition 
     265// the user requests a polynomial of order Norder
     266psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
     267                                    psPolynomialType type)
     268{
     269    PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
     270    PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
     271
     272    psS32 x = 0;
     273    psS32 y = 0;
     274    psS32 nXterm = nXorder + 1;
     275    psS32 nYterm = nYorder + 1;
     276    psPolynomial2D* newPoly = NULL;
     277
     278    newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
     279    // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
     280    // XXX EAM : me, being lazy
     281
     282    newPoly->type = type;
     283    newPoly->nX = nXorder;
     284    newPoly->nY = nYorder;
     285
     286    newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
     287    newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
     288    newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
     289    for (x = 0; x < nXterm; x++) {
     290        newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
     291        newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
     292        newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
     293    }
     294    for (x = 0; x < nXterm; x++) {
     295        for (y = 0; y < nYterm; y++) {
     296            newPoly->coeff[x][y] = 0.0;
     297            newPoly->coeffErr[x][y] = 0.0;
     298            newPoly->mask[x][y] = 0;
     299        }
     300    }
     301    return(newPoly);
     302}
     303
     304// XXX EAM : BuildSums2D in analogy with BuildSums1D
     305static psImage *BuildSums2D(psImage* sums,
     306                            psF64 x,      psF64 y,
     307                            psS32 nXterm, psS32 nYterm)
     308{
     309    psS32 nXsum = 0;
     310    psS32 nYsum = 0;
     311    psF64 xSum = 1.0;
     312    psF64 ySum = 1.0;
     313
     314    nXsum = 2*nXterm;
     315    nYsum = 2*nYterm;
     316    if (sums == NULL) {
     317        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
     318    }
     319    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
     320        psFree (sums);
     321        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
     322    }
     323
     324    ySum = 1.0;
     325    for (int j = 0; j < nYsum; j++) {
     326        xSum = ySum;
     327        for (int i = 0; i < nXsum; i++) {
     328            sums->data.F64[i][j] = xSum;
     329            xSum *= x;
     330        }
     331        ySum *= y;
     332    }
     333    return (sums);
     334}
     335
    293336// XXX EAM : test version of 2d fitting
    294337psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
    … …  
    387430}   
    388431
    389 // write out the terms of the given 1D polynomial
    390 void psPolynomial1DDump (psPolynomial1D *poly) {
    391 
    392     for (int i = 0; i < poly->n + 1; i++) {
    393         fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]);
    394     }
    395 }   
    396 
    397432psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
    398433                            const psVector* x,

  • trunk/psphot/src/psUtils.c

    r4251 r4375  
    448448}
    449449
     450psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) {
     451
     452  int nBin = (upper - lower) / delta;
     453
     454  psVector *out = psVectorAlloc (nBin, type);
     455 
     456  for (int i = 0; i < nBin; i++) {
     457    out->data.F64[i] = lower + i * delta;
     458  }
     459
     460  return (out);
     461}

  • trunk/psphot/src/psphot-utils.c

    r4251 r4375  
    7777        params = model->params;
    7878        dparams = model->dparams;
    79         fprintf (f, "%7.1f %7.1f  %5.1f %7.1f  %7.4f  ",
     79        fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f  ",
    8080                 params[0].data.F32[2], params[0].data.F32[3],
    81                  params[0].data.F32[0], params[0].data.F32[1], (dparams[0].data.F32[1]/params[0].data.F32[1]));
    82         for (j = 0; j < model->params->n - 4; j++) {
    83             fprintf (f, "%9.3g ", params[0].data.F32[j+4]);
    84         }
    85         for (j = 0; j < model->params->n - 4; j++) {
    86             fprintf (f, "%9.3g ", dparams[0].data.F32[j+4]);
     81                 params[0].data.F32[0], -2.5*log10(params[0].data.F32[1]), (dparams[0].data.F32[1]/params[0].data.F32[1]));
     82        for (j = 0; j < model->params->n - 4; j++) {
     83            fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
     84        }
     85        for (j = 0; j < model->params->n - 4; j++) {
     86            fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
    8787        }
    8888        fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
     89    }
     90    fclose (f);
     91    return true;
     92}
     93
     94// dump the sources to an output file
     95bool DumpModelFLT (psArray *sources, char *filename)
     96{
     97
     98    int i, j;
     99    FILE *f;
     100    psVector *params;
     101    psVector *dparams;
     102    psModel  *model;
     103
     104    f = fopen (filename, "w");
     105    if (f == NULL) {
     106        psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
     107        return false;
     108    }
     109
     110    // write sources with models first
     111    for (i = 0; i < sources->n; i++) {
     112        psSource *source = (psSource *) sources->data[i];
     113        model = (psModel  *) source->modelFLT;
     114        if (model == NULL) continue;
     115        params = model->params;
     116        dparams = model->dparams;
     117        fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f  ",
     118                 params[0].data.F32[2], params[0].data.F32[3],
     119                 params[0].data.F32[0], -2.5*log10(params[0].data.F32[1]), (dparams[0].data.F32[1]/params[0].data.F32[1]));
     120        for (j = 0; j < model->params->n - 4; j++) {
     121            fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
     122        }
     123        for (j = 0; j < model->params->n - 4; j++) {
     124            fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
     125        }
     126        fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
     127    }
     128    fclose (f);
     129    return true;
     130}
     131
     132// dump the sources to an output file
     133bool DumpModelNULL (psArray *sources, char *filename)
     134{
     135
     136    int i, j;
     137    FILE *f;
     138    psVector *params;
     139    psVector *dparams;
     140    psModel  *model;
     141
     142    f = fopen (filename, "w");
     143    if (f == NULL) {
     144        psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
     145        return false;
     146    }
     147
     148    // write sources with models first
     149    for (i = 0; i < sources->n; i++) {
     150        psSource *source = (psSource *) sources->data[i];
     151        if (source->modelFLT != NULL) continue;
     152        if (source->modelPSF != NULL) continue;
     153        if (source->moments == NULL) continue;
     154        fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n",
     155                 source->peak->x, source->peak->y, source->peak->counts,
     156                 source->moments->x, source->moments->y,
     157                 source->moments->Sx, source->moments->Sy,
     158                 source->moments->Sum, source->moments->Peak,
     159                 source->moments->Sky, source->moments->SN,
     160                 source->moments->nPixels, source->type);
    89161    }
    90162    fclose (f);
    … …  
    118190                 params[0].data.F32[0], params[0].data.F32[1]);
    119191        for (j = 0; j < model->params->n - 4; j++) {
    120             fprintf (f, "%9.3g ", params[0].data.F32[j+4]);
    121         }
    122         for (j = 0; j < model->params->n - 4; j++) {
    123             fprintf (f, "%9.3g ", dparams[0].data.F32[j+4]);
     192            fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
     193        }
     194        for (j = 0; j < model->params->n - 4; j++) {
     195            fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
    124196        }
    125197        fprintf (f, ":  %7.1f (%d %d)\n", model[0].chisq, model[0].nDOF, model[0].nIter);
    … …  
    138210}
    139211
     212bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask) {
     213
     214    float obsSum = 0;
     215    float fitSum = 0;
     216    float sky = model->params->data.F32[0];
     217
     218    psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type);
     219    fitSum = modelFluxFunc (model->params);
     220
     221    for (int ix = 0; ix < image->numCols; ix++) {
     222        for (int iy = 0; iy < image->numRows; iy++) {
     223            if (mask->data.U8[iy][ix]) continue;
     224            obsSum += image->data.F32[iy][ix] - sky;
     225            // fitSum += psModelEval (model, image, ix, iy) - sky;
     226        }
     227    }
     228    if (obsSum <= 0) return false;
     229    if (fitSum <= 0) return false;
     230
     231    *fitMag = -2.5*log10(fitSum);
     232    *obsMag = -2.5*log10(obsSum);
     233    return (true);
     234}
     235
    140236# if (0)
    141237void pmSourceMaskSaturated (psSource *source, float saturate) {

  • trunk/psphot/src/psphot.c

    r4251 r4375  
    2929
    3030    // use stellar objects SN > PSF_SN_LIM
    31     psf = choose_psf_model (config, sources);
     31    psf = choose_psf_model (config, sources, sky);
    3232
    3333    // identify defects based on 1-pixel wide x and y direction 2nd moments
    34     find_defects (zap, sources, config, psf);
     34    // find_defects (zap, sources, config, psf);
    3535
    3636    // test PSF on all except SATURATE and DEFECT (artifacts)
    … …  
    4646 
    4747    DumpImage (imdata->image, argv[2]);
    48     DumpModelPSF (sources, "sources.dat");
     48    DumpModelPSF (sources, "psfsources.dat");
     49    DumpModelFLT (sources, "fltsources.dat");
     50    DumpModelNULL (sources, "nullsources.dat");
     51    DumpMoments (sources, "moments.dat");
     52
    4953    exit (0);
    5054}

  • trunk/psphot/src/psphot.h

    r4251 r4375  
    55# include <unistd.h>
    66# include <stdlib.h>
     7# include <math.h>
     8
     9# define  M_PI 3.14159265358979323846264338328      /* pi */
    710
    811typedef struct {
    … …  
    2225  psVector   *mask;
    2326  psVector   *metric;
    24   psStats    *metricStats;
     27  psVector   *fitMag;
     28  float       ApResid;
     29  float       dApResid;
     30  float       skyBias;
    2531} pmPSF_Test;
    2632
    27 bool psTimerStart (char *name);
    28 psF64 psTimerMark (char *name);
     33// used by mask operations
     34# define AND true
     35# define OR false
    2936
    30 bool psImageInit (psImage *image,...);
    31 void psImageSmooth (psImage *image, float sigma, float Nsigma);
    32 psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
    3337
     38// top-level psphot functions
     39psImageData *setup (psMetadata *config);
     40psStats *image_stats (psImageData *imdata, psMetadata *config);
     41psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky);
     42psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks);
     43pmPSF *choose_psf_model (psMetadata *config, psArray *sources, psStats *sky);
     44bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
     45bool test_psf_scatter (psArray *sources);
     46bool mark_psf_sources (psArray *sources, psMetadata *config);
     47bool subtract_psf_sources (psArray *sources);
     48bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
     49bool subtract_galaxies (psArray *sources, psMetadata *config);
     50bool subtract_psf_source (psSource *source);
     51bool mark_psf_source (psSource *source, float shapeNsigma);
     52bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
     53
     54// psf utilities
     55pmPSF *pmPSFAlloc (psModelType type);
     56pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName);
     57pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
     58bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
     59psModel *psModelFromPSF (psModel *model, pmPSF *psf);
     60bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask);
     61bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS);
     62
     63// output options
    3464bool DumpPeaks (psArray *sources, char *filename);
    3565bool DumpMoments (psArray *sources, char *filename);
    … …  
    3868bool DumpImage (psImage *image, char *filename);
    3969
     70// psphot utilities
    4071int usage ();
    41 int get_argument (int argc, char **argv, char *arg);
    42 int remove_argument (int N, int *argc, char **argv);
    4372psMetadata *load_args (int *argc, char **argv);
    44 
     73int by_SN (const void **a, const void **b);
    4574psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid);
    4675psArray *SelectPSFStars (psArray *sources);
    47 
    48 psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
    49 
    50 psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
    51 psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
    52 psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
    53 psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
    54 psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
    55 void psPolynomial2DDump (psPolynomial2D *poly);
    56 psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
    57 psPolynomial2D* RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
    58 psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
    59 psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
    60 void psPolynomial1DDump (psPolynomial1D *poly);
    61 
    62 psS32 psLogArguments (int *argc, char **argv);
    63 psS32 psTraceArguments (int *argc, char **argv);
    64 psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row);
    65 psVector *psGetRowVectorFromImage(psImage *image, psU32 row);
    66 
    67 pmPSF *pmPSFAlloc (psModelType type);
    68 pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName);
    69 pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
    70 
    71 bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
    72 psModel *psModelFromPSF (psModel *model, pmPSF *psf);
    73 void pmSourceMaskRegion (psSource *source, psRegion *region);
    74 void pmSourceMaskSaturated (psSource *source, float saturate);
    75 
    76 // top-level psphot functions
    77 psImageData *setup (psMetadata *config);
    78 psStats *image_stats (psImageData *imdata, psMetadata *config);
    79 psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky);
    80 psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks);
    81 pmPSF *choose_psf_model (psMetadata *config, psArray *sources);
    82 bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
    83 bool test_psf_scatter (psArray *sources);
    84 bool mark_psf_sources (psArray *sources, psMetadata *config);
    85 bool subtract_psf_sources (psArray *sources);
    86 bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
    87 bool subtract_galaxies (psArray *sources, psMetadata *config);
    88 
    89 int by_SN (const void **a, const void **b);
    90 bool subtract_psf_source (psSource *source);
    91 bool mark_psf_source (psSource *source, float shapeNsigma);
    92 bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
    93 
    9476bool pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
    9577bool pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius);
    96 psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius);
     78
     79// fitsource utilities
     80bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y);
     81psMetadata *fs_args (int *argc, char **argv);
     82int fs_usage ();
    9783
    9884// image mask functions
    … …  
    10288void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue);
    10389
    104 # define AND true
    105 # define OR false
     90// basic image functions
     91bool psImageInit (psImage *image,...);
     92void psImageSmooth (psImage *image, float sigma, float Nsigma);
     93psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
     94psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius);
     95
     96// misc global functions
     97psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row);
     98psVector *psGetRowVectorFromImage(psImage *image, psU32 row);
     99
     100// utilities
     101bool psTimerStart (char *name);
     102psF64 psTimerMark (char *name);
     103psS32 psLogArguments (int *argc, char **argv);
     104psS32 psTraceArguments (int *argc, char **argv);
     105int get_argument (int argc, char **argv, char *arg);
     106int remove_argument (int N, int *argc, char **argv);
     107psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
     108psVector *psVectorCreate (double lower, double upper, double delta, psElemType type);
     109
     110// polynomial functions
     111psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
     112psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
     113psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
     114psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
     115void psPolynomial2DDump (psPolynomial2D *poly);
     116psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
     117psPolynomial2D* RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
     118psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
     119
     120psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
     121void psPolynomial1DDump (psPolynomial1D *poly);
     122psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly);
     123psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
     124psPolynomial1D* Polynomial1DAlloc(psS32 nOrder, psPolynomialType type);
     125psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
     126
     127// deprecated
     128void pmSourceMaskRegion (psSource *source, psRegion *region);
     129void pmSourceMaskSaturated (psSource *source, float saturate);

  • trunk/psphot/src/pspsf.c

    r4251 r4375  
    3434    test->modelFLT    = psArrayAlloc (sources->n);
    3535    test->modelPSF    = psArrayAlloc (sources->n);
    36     test->metricStats = NULL;
    3736    test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
     37    test->fitMag      = psVectorAlloc (sources->n, PS_TYPE_F64);
    3838    test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
     39    test->ApResid     = 0;
     40    test->dApResid    = 0;
     41    test->skyBias     = 0;
    3942
    4043    for (int i = 0; i < test->modelFLT->n; i++) {
     44        test->mask->data.U8[i]  = 0;
    4145        test->modelFLT->data[i] = NULL;
    4246        test->modelPSF->data[i] = NULL;
    43         test->mask->data.U8[i]  = 0;
     47        test->metric->data.F64[i] = 0;
     48        test->fitMag->data.F64[i] = 0;
    4449    }   
    4550    return (test);
    4651}
    4752
     53// test->mask values indicate reason source was rejected:
     54// 1: outlier in psf polynomial fit
     55// 2: flt model failed to converge
     56// 3: psf model failed to converge
     57// 4: invalid source photometry
    4858pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS)
    4959{
    50  
     60    bool status;
     61    float obsMag;
     62    float fitMag;
    5163    float x;
    5264    float y;
    53     bool status;
    5465    int Nflt = 0;
    5566    int Npsf = 0;
    56     psVector *metricMask = psVectorAlloc (sources->n, PS_TYPE_U8);
     67
    5768    pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
    5869
    … …  
    6677
    6778        // set temporary object mask and fit object
    68         // fit model as FLT, not PSF (mask & skip poor fits)
     79        // fit model as FLT, not PSF
    6980        psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
    7081        status = pmSourceFitModel (source, model, false);
    7182        psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
    7283
     84        // exclude the poor fits
    7385        if (!status) {
    74           test->mask->data.U8[i] = 1;
     86          test->mask->data.U8[i] = 2;
    7587          continue;
    7688        }
    … …  
    8496    // stage 2: construct a psf (pmPSF) from this collection of model fits
    8597    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
    86    
    87     // count valid sources
    88     int Nkeep = 0;
    89     for (int i = 0; i < sources->n; i++) {
    90       if (test->mask->data.U8[i]) continue;
    91       Nkeep++;
    92     }
    93     psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates\n", Nkeep, sources->n);
    9498
    9599    // stage 3: refit with fixed shape parameters
    96100    psTimerStart ("fit");
    97101    for (int i = 0; i < test->sources->n; i++) {
     102        // masked for: bad model fit, outlier in parameters
     103        if (test->mask->data.U8[i]) continue;
     104
    98105        psSource *source = test->sources->data[i];
    99106        psModel  *modelFLT = test->modelFLT->data[i];
    100107
    101         // masked for: bad model fit, outlier in parameters
    102         if (test->mask->data.U8[i]) continue;
    103         if (modelFLT == NULL) {
    104           psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
    105           continue;
    106         }
    107         psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
     108        // set shape for this model based on PSF
     109        psModel *modelPSF = psModelFromPSF (modelFLT, test->psf);
    108110        x = source->peak->x;
    109111        y = source->peak->y;
    … …  
    111113        psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
    112114        status = pmSourceFitModel (source, modelPSF, true);
    113         psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
    114115
    115116        // skip poor fits
    116117        if (!status) {
    117           test->mask->data.U8[i] = 1;
    118           continue;
     118            test->mask->data.U8[i] = 3;
     119            goto next_source;
    119120        }
     121
     122        // otherwise, save the resulting model
    120123        test->modelPSF->data[i] = modelPSF;
     124
     125        // XXX : use a different aperture radius from the fit radius?
     126        // XXX : use a different estimator for the local sky?
     127        // XXX : pass 'source' as input?
     128        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
     129            test->mask->data.U8[i] = 4;
     130            goto next_source;
     131        }           
     132
     133        test->metric->data.F64[i] = obsMag - fitMag;
     134        test->fitMag->data.F64[i] = fitMag;
    121135        Npsf ++;
     136
     137    next_source:
     138        psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
     139
    122140    }
    123141    psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    … …  
    125143    DumpModelFits (test->modelPSF, "modelsPSF.dat");
    126144
    127     // stage 4: measure PSF metric (Ap-Fit)
    128     // should I be saving the aperture photometry for all sources?
    129     // save the metric for all stars
    130     for (int i = 0; i < test->sources->n; i++) {
    131 
    132         float obsSum = 0;
    133         float fitSum = 0;
    134      
    135         // is this metricMask redundant with test->mask ?
    136         metricMask->data.U8[i] = 1;
    137         test->metric->data.F64[i] = 0;
    138 
    139         if (test->mask->data.U8[i]) continue;
    140 
    141         psModel *model = test->modelPSF->data[i];
    142         if (model == NULL) {
    143           psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
    144           continue;
    145         }
    146 
    147         psImage *image = ((psSource *)test->sources->data[i])->pixels;
    148         psImage *mask  = ((psSource *)test->sources->data[i])->mask;
    149 
    150         // this metric is Ap-Fit
    151         float sky = model->params->data.F32[0];
    152         for (int ix = 0; ix < image->numCols; ix++) {
    153             for (int iy = 0; iy < image->numRows; iy++) {
    154                 if (mask->data.U8[iy][ix]) continue;
    155                 obsSum += image->data.F32[iy][ix] - sky;
    156                 fitSum += psModelEval (model, image, ix, iy) - sky;
    157             }
    158         }
    159         // fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
    160         // keep the good metrics
    161         if ((fitSum > 0) && (obsSum > 0)) {
    162             metricMask->data.U8[i] = 0;
    163             test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
    164         } else {
    165           test->mask->data.U8[i] = 1;
    166         }
    167     }
    168 
    169     // XXX use robust stats, not clipped stats
    170     psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    171     psVectorStats (stats, test->metric, NULL, metricMask, 1);
    172     test->metricStats = stats;
    173     psLogMsg ("psphot.pspsf", 3, "test model %s, metric: %f +/- %f\n", modelName, stats->clippedMean, stats->clippedStdev);
     145    // XXX this function wants aperture radius from pmSourcePhotometry
     146    pmPSFMetricModel (test, RADIUS);
     147    psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n",
     148              modelName, test->ApResid, test->dApResid, test->skyBias);
    174149    return (test);
    175150}
    … …  
    181156    // construct the fit vectors from the collection of objects
    182157    // use the mask to ignore missing fits
    183     psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
    184     psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
    185     psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
     158    psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
     159    psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
     160    psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
    186161    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
    187162
    … …  
    208183            // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
    209184        }
    210         // if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
    211185
    212186        psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
     187        // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
    213188        // psPolynomial2DDump (psf->params->data[i]);
    214 
    215         // count valid sources
    216         int Nkeep = 0;
    217         for (int j = 0; j < mask->n; j++) {
    218           if (mask->data.U8[j]) continue;
    219           Nkeep++;
    220         }
    221         psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
    222 
    223189    }
    224190    return (true);
    … …  
    241207}
    242208
     209bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
     210
     211  float dBin;
     212  int   nKeep, nSkip;
     213 
     214
     215  // the measured (aperture - fit) magnitudes (dA == test->metric)
     216  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
     217  //     dA = dAo + dsky/flux
     218  //   where flux is the flux of the star
     219  // we fit this trend to find the infinite flux aperture correction (dAo),
     220  //   the nominal sky bias (dsky), and the error on dAo
     221  // the values of dA are contaminated by stars with close neighbors in the aperture
     222  //   we use an outlier rejection to avoid this bias
     223
     224  // rflux = ten(0.4*fitMag);
     225  psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
     226  for (int i = 0; i < test->sources->n; i++) {
     227    if (test->mask->data.U8[i]) continue;
     228    rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
     229  }
     230  // psScalar *t1 = psScalar(0.4);
     231  // psVector *v1 = psBinaryOp (NULL, test->fitMag, "*", t1);
     232  // psVector *rflux = psUnaryOp (NULL, v1, "ten");
     233  // psFree (t1);
     234  // psFree (v1);
     235
     236  // find min and max of (1/flux):
     237  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
     238  psVectorStats (stats, rflux, NULL, test->mask, 0xff);
     239 
     240  // build binned versions of rflux, metric
     241  dBin = (stats->max - stats->min) / 10.0;
     242  psVector *rfBin = psVectorCreate (stats->min, stats->max, dBin, PS_TYPE_F64);
     243  psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
     244  psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
     245  psFree (stats);
     246
     247  psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
     248
     249  for (int i = 0; i < daBin->n; i++) {
     250
     251    psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
     252    tmp->n = 0;
     253
     254    // accumulate data within bin range
     255    for (int j = 0; j < test->sources->n; j++) {
     256      // masked for: bad model fit, outlier in parameters
     257      if (test->mask->data.U8[j]) continue;
     258   
     259      // skip points with extreme dA values
     260      if (fabs(test->metric->data.F64[j]) > 0.5) continue;
     261
     262      // skip points outside of this bin
     263      if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue;
     264      if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
     265
     266      tmp->data.F64[tmp->n] = test->metric->data.F64[j];
     267      tmp->n ++;
     268    }
     269
     270    // is this a valid point?
     271    maskB->data.U8[i] = 0;
     272    if (tmp->n < 2) {
     273      maskB->data.U8[i] = 1;
     274      continue;
     275    }
     276
     277    // dA values are contaminated with low outliers
     278    // measure statistics only on upper 50% of points
     279    psVectorSort (tmp, tmp);
     280    nKeep = 0.5*tmp->n;
     281    nSkip = tmp->n - nKeep;
     282
     283    psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
     284    for (int j = 0; j < tmp2->n; j++) {
     285      tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
     286    }
     287
     288    stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     289    psVectorStats (stats, tmp2, NULL, NULL, 0);
     290    psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
     291
     292    daBin->data.F64[i] = stats->sampleMedian;
     293  }
     294
     295  // linear fit to rfBin, daBin
     296  psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
     297  poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
     298
     299  psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
     300  psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
     301
     302  stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
     303  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
     304
     305  test->ApResid = poly->coeff[0];
     306  test->dApResid = stats->clippedStdev;
     307  test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     308
     309  return true;
     310}

  • trunk/psphot/src/setup.c

    r4251 r4375  
    1818    psTimerStart ("psphot");
    1919
    20     // setup header template, specify COMMENT and HISTORY as multis
     20    // setup header template, specify COMMENT and HISTORY as MULTI
    2121    header = psMetadataAlloc ();
    2222    psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
    … …  
    8989        }
    9090    }
     91    DumpImage (mask, "mask.fits");
    9192
    9293    psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));

  • trunk/psphot/src/source_moments.c

    r4251 r4375  
    55    bool     status  = false;
    66    psArray *sources = NULL;
     7    pmPSFClump psfClump;
    78
    89    psTimerStart ("psphot");
    … …  
    2526
    2627        // this should use ROBUST not SAMPLE median, but it is broken
    27         pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
    28         pmSourceMoments (source, RADIUS);
     28        status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     29        if (!status) continue;
     30
     31        status = pmSourceMoments (source, RADIUS);
     32        if (!status) continue;
     33       
    2934        psArrayAdd (sources, 100, source);
    3035    }
    3136
    3237    // group into STAR, COSMIC, GALAXY, SATURATED
    33     pmSourceRoughClass (sources, config);
     38    psfClump = pmSourcePSFClump (sources, config);
     39    pmSourceRoughClass (sources, config, psfClump);
    3440
    3541    // make this optional
    36     DumpMoments (sources, "moments.dat");
     42    // DumpMoments (sources, "moments.dat");
    3743
    3844    psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot"));
    … …  
    4147    return (sources);
    4248}
     49
     50// filter out bad peaks (eg, no valid pixels available for sky)
Note: See TracChangeset for help on using the changeset viewer.