IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15000 for trunk/psphot


Ignore:
Timestamp:
Sep 24, 2007, 11:27:58 AM (19 years ago)
Author:
eugene
Message:

update from eam_branch_20070921

Location:
trunk/psphot
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/doc/psfmodel.txt

    r9772 r15000  
     1
     22007.09.21
     3
     4  there are three places where we can choose to use errors in the fits or not:
     5
     6  * non-linear fitting of the models to the pixel flux distribution (poissonErrorsPhotLMM)
     7  * linear fitting of the models to the pixel flux distribution (poissonErrorsPhotLin)
     8  * fitting of the 2D variations in the psf parameters (poissonErrorsParams)
     9  * fitting of the 2D variations in the aperture residuals
     10
     112007.09.20
     12
     13  I am upgrading the PSF model to allow the parameter variation to be
     14  modeled with pmTrend2D instead of just polynomials.  I am making a
     15  list of places to modify the code:
     16
     17pmPSFAlloc : need a method beyong psfTrendMask to carry in the psf
     18options
     19
     20pmPSF_ModelToFit : no need to change these
     21
     22update pmPSFBuildSimple to set the parameters of the pmTrend, which
     23ever is used.
     24
     25pmPSFtry.c: some significant re-work!
     26
     27
     28
     29pmPSF_IO : need new functions to save / load the trend (psImages)
    130
    2312006.10.27

  • trunk/psphot/src/models/pmModel_STRAIL.c

    r14655 r15000  
    530530 
    531531    for (int i = 4; i < 7; i++) {
    532       psPolynomial2D *poly = psf->params->data[i-4];
    533         out[i] = psPolynomial2DEval (poly, out[2], out[3]);
     532      pmTrend2D *trend = psf->params->data[i-4];
     533        out[i] = pmTrend2DEval (trend, out[2], out[3]);
    534534    }
    535535    return(true);
    … …  
    554554    for (int i = 0; i < psf->params->n; i++) {
    555555        if (i == PM_PAR_SKY) continue;
    556         psPolynomial2D *poly = psf->params->data[i];
    557         assert (poly);
    558         PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
     556        pmTrend2D *trend = psf->params->data[i];
     557        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
    559558    }
    560559

  • trunk/psphot/src/models/pmModel_TEST1.c

    r14655 r15000  
    217217            out[i] = in[i];
    218218        } else {           
    219             psPolynomial2D *poly = psf->params->data[i];
    220             out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
     219            pmTrend2D *trend = psf->params->data[i];
     220            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
    221221        }
    222222    }
    … …  
    246246    for (int i = 0; i < psf->params->n; i++) {
    247247        if (i == PM_PAR_SKY) continue;
    248         psPolynomial2D *poly = psf->params->data[i];
    249         assert (poly);
    250         PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
     248        pmTrend2D *trend = psf->params->data[i];
     249        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
    251250    }
    252251

  • trunk/psphot/src/psphot.h

    r14963 r15000  
    4747bool            psphotReplaceAll (psArray *sources, psMaskType maskVal);
    4848bool            psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark);
    49 bool            psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup);
     49bool            psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup);
    5050bool            psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag);
    5151bool            psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background, psMaskType maskVal, psMaskType mark);

  • trunk/psphot/src/psphotApResid.c

    r14965 r15000  
    11# include "psphotInternal.h"
    22
     3# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
    34// measure the aperture residual statistics and 2D variations
    45bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark)
    … …  
    6364        model = source->modelPSF;
    6465
    65         if (source->type != PM_SOURCE_TYPE_STAR) continue;
    66         if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
    67         if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
    68         if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
    69         if (source->mode &  PM_SOURCE_MODE_POOR) continue;
     66        if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
     67        if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
     68        if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
     69        if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
     70        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
    7071
    7172        // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
    … …  
    7374        if (!pmSourceMagnitudes (source, psf, photMode, maskVal, mark)) {
    7475            Nskip ++;
     76            psTrace ("psphot", 3, "skip : bad source mag");
    7577            continue;
    7678        }
    7779
    78         if (!isfinite(source->apMag) || !isfinite(source->apMag)) {
     80        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    7981            Nfail ++;
     82            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
    8083            continue;
    8184        }
    … …  
    8790        if ((MAX_AP_OFFSET > 0) && (fabs(dap) > MAX_AP_OFFSET)) {
    8891            Nfail ++;
     92            psTrace ("psphot", 3, "fail : bad dap %f %f", dap, MAX_AP_OFFSET);
    8993            continue;
    9094        }
    … …  
    99103        dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    100104
     105        psVectorExtend (mag,     100, 1);
    101106        psVectorExtend (mask,    100, 1);
    102107        psVectorExtend (xPos,    100, 1);
    … …  
    146151    }
    147152
     153    // XXX catch error condition
    148154    psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
    149155
    … …  
    152158    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
    153159    psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
    154 
    155     psphotSaveImage (NULL, psf->ApTrend->map->map, "apresid.map.fits");
    156160
    157161    if (psTraceGetLevel("psphot") >= 5) {
    … …  
    170174    float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
    171175    float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
     176
    172177    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
    173178    psf->dApResid = errorFloor;
    … …  
    206211*/
    207212
    208 bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup) {
     213bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
    209214
    210215    psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    … …  
    225230    psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    226231    psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
     232    psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_U8);
    227233
    228234    int n = 0;
    … …  
    233239            dMSubset->data.F32[j] = dMag->data.F32[N];
    234240            dASubset->data.F32[j] = dap->data.F32[N];
     241            mkSubset->data.U8[j]  = mask->data.U8[N];
    235242        }
    236243        dMSubset->n = j;
    237244        dASubset->n = j;
     245        mkSubset->n = j;
    238246
    239247        psStatsInit (statsS);
    240248        psStatsInit (statsM);
    241249
    242         psVectorStats (statsS, dASubset, NULL, NULL, 0xff);
    243         psVectorStats (statsM, dMSubset, NULL, NULL, 0xff);
     250        psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff);
     251        psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff);
    244252
    245253        dSo->data.F32[i] = statsS->robustStdev;
    … …  
    250258    psFree (dMSubset);
    251259    psFree (dASubset);
     260    psFree (mkSubset);
    252261   
    253262    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    … …  
    286295    }
    287296
     297    // the mask marks the values not used to calculate the ApTrend
     298    psVectorInit (mask, 0);
     299
    288300    // XXX stats structure for use by ApTrend : make parameters user setable
    289301    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    … …  
    305317    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
    306318    int nGroup = PS_MAX (3*Nx*Ny, 10);
    307     psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, nGroup);
     319    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
    308320
    309321    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);

  • trunk/psphot/src/psphotChoosePSF.c

    r14966 r15000  
    1111    sources = psArraySort (sources, pmSourceSortBySN);
    1212
     13    // structure to store user options defining the psf
     14    pmPSFOptions *options = pmPSFOptionsAlloc ();
     15
     16    // load user options from the recipe. no need to check existence -- they are
    1317    // array to store candidate PSF stars
    1418    int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS");
    15     PS_ASSERT (status, NULL);
     19    assert (status);
    1620
    1721    float PSF_MIN_DS = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS");
    18     PS_ASSERT (status, NULL);
     22    assert (status);
    1923
    2024    // supply the measured sky variance for optional constant errors (non-poissonian)
    2125    float SKY_SIG = psMetadataLookupF32 (&status, recipe, "SKY_SIG");
    22     PS_ASSERT (status, NULL);
     26    assert (status);
    2327
    2428    // use poissonian errors or local-sky errors
    25     bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
    26     PS_ASSERT (status, NULL);
     29    options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LMM");
     30    assert (status);
    2731
    2832    // use poissonian errors or local-sky errors
    29     bool PSF_PARAM_WEIGHTS = psMetadataLookupBool (&status, recipe, "PSF_PARAM_WEIGHTS");
    30     PS_ASSERT (status, NULL);
     33    options->poissonErrorsPhotLin = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LIN");
     34    assert (status);
     35
     36    // use poissonian errors or local-sky errors
     37    options->poissonErrorsParams = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PARAMS");
     38    assert (status);
    3139
    3240    // how to model the PSF variations across the field
    33     psMetadata *md = psMetadataLookupMetadata (&status, recipe, "PSF.TREND.MASK");
    34     PS_ASSERT (status, NULL);
     41    options->psfTrendMode = pmTrend2DModeFromString (psMetadataLookupStr (&status, recipe, "PSF.TREND.MODE"));
     42    assert (status);
     43   
     44    options->psfTrendNx = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NX");
     45    assert (status);
     46
     47    options->psfTrendNy = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NY");
     48    assert (status);
    3549
    3650    // get the fixed PSF fit radius
    3751    // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
    38     float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    39     PS_ASSERT (status, NULL);
    40 
    41     pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), POISSON_ERRORS);
    42 
    43     psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md);
    44     if (!psfTrendMask) {
    45         psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
    46         return NULL;
    47     }
     52    options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
     53    assert (status);
     54
     55    options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     56
     57    // dimensions of the field for which the PSF is defined
     58    options->psfFieldNx = readout->image->numCols;
     59    options->psfFieldNy = readout->image->numRows;
     60    options->psfFieldXo = readout->image->col0;
     61    options->psfFieldYo = readout->image->row0;
     62
     63    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM);
    4864
    4965    psArray *stars = psArrayAllocEmpty (sources->n);
    … …  
    6682        psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates");
    6783        psFree (stars);
    68         psFree (psfTrendMask);
    6984        return NULL;
    7085    }
    … …  
    91106        psMetadataItem *item = psListGetAndIncrement (iter);
    92107        char *modelName = item->data.V;
    93         models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark);
     108        models->data[i] = pmPSFtryModel (stars, modelName, options, maskVal, mark);
    94109    }
    95110
    … …  
    97112    psFree (list);
    98113    psFree (stars);
    99     psFree (psfTrendMask);
    100114
    101115    // select the best of the models
    … …  
    128142    if (psTraceGetLevel("psphot") >= 5) {
    129143        for (int i = PM_PAR_SXX; i < try->psf->params->n; i++) {
    130             psPolynomial2D *poly = try->psf->params->data[i];
    131             for (int nx = 0; nx <= poly->nX; nx++) {
    132                 for (int ny = 0; ny <= poly->nY; ny++) {
    133                     if (poly->mask[nx][ny]) continue;
    134                     fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny);
    135                 }
    136             }
     144            // XXX re-write the output or some other status info
    137145        }
    138146    }
    … …  
    236244            // use pmModelSub because modelFlux has not been generated
    237245            assert (source->maskObj);
    238             psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
     246            psImageKeepCircle (source->maskObj, x, y, options->radius, "OR", PM_MASK_MARK);
    239247            pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
    240             psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
     248            psImageKeepCircle (source->maskObj, x, y, options->radius, "AND", PS_NOT_U8(PM_MASK_MARK));
    241249        }
    242250
    … …  
    283291        pmSourcesWritePSFs (try->sources, "psfstars.dat");
    284292        pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
    285         psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
    286         psMetadataConfigWrite (psfData, "psfmodel.dat");
    287         psFree (psfData);
     293        // XXX need alternative output function
     294        // psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
     295        // psMetadataConfigWrite (psfData, "psfmodel.dat");
    288296        psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
    289297        exit (0);
    … …  
    306314    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
    307315
     316    psFree (options);
    308317    return (psf);
    309318}

  • trunk/psphot/src/psphotRoughClass.c

    r14760 r15000  
    11# include "psphotInternal.h"
     2# define CHECK_STATUS(S,MSG) { \
     3  if (!status) { \
     4    psError(PSPHOT_ERR_CONFIG, false, "missing PSF Clump entry: %s\n", MSG); \
     5    return false; \
     6  } }
    27
    38// 2006.02.02 : no leaks
    49bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool havePSF, psMaskType maskSat) {
    510
    6     static pmPSFClump   psfClump;
     11    bool status;
     12    static pmPSFClump psfClump;
    713
    814    psTimerStart ("psphot");
    … …  
    1117        // determine the PSF parameters from the source moment values
    1218        psfClump = pmSourcePSFClump (sources, recipe);
     19        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X", 0, "psf clump center", psfClump.X);
     20        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y", 0, "psf clump center", psfClump.Y);
     21        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", 0, "psf clump center", psfClump.dX);
     22        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", 0, "psf clump center", psfClump.dY);
    1323    } else {
    1424        // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
    15         bool status_x, status_y, status_a;
    16         float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
    17         float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
    18         float ANGLE  = psMetadataLookupF32 (&status_a, recipe, "ANGLE");
    19         if (!status_x | !status_y | !status_a) {
    20             psError(PSPHOT_ERR_CONFIG, false, "FWHM_X, FWHM_Y, or ANGLE not defined");
    21             return false;
    22         }
    23 
    24         psEllipseAxes axes;
    25         axes.major = FWHM_X / (2.0*sqrt(2.0*log(2.0)));
    26         axes.minor = FWHM_Y / (2.0*sqrt(2.0*log(2.0)));
    27         axes.theta = ANGLE;
    28         psEllipseShape shape = psEllipseAxesToShape (axes);
    29 
    30         psfClump.X   = shape.sx;
    31         psfClump.Y   = shape.sy;
    32         psfClump.dX  = 0.1*shape.sx;
    33         psfClump.dY  = 0.1*shape.sy;
    34         // dX,dY are somewhat crudely defined, but only used to select PSF candidates.
    35         // if we already have a PSF, this is not actually used...
     25        psfClump.X  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");  CHECK_STATUS (status, "PSF.CLUMP.X");
     26        psfClump.Y  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y");  CHECK_STATUS (status, "PSF.CLUMP.Y");
     27        psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX");  CHECK_STATUS (status, "PSF.CLUMP.DX");
     28        psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");  CHECK_STATUS (status, "PSF.CLUMP.DY");
    3629    }
    3730

  • trunk/psphot/src/psphotSourcePlots.c

    r13900 r15000  
    2020    int dY = 0;                         // height of row so far
    2121    int NX = 20*DX;                     // full width of output image
    22     int NY = 0;                         // total height of output image
     22    int NY = DY;                        // total height of output image so far
    2323
    24     // first, examine the PSF and SAT stars:
    25     // - determine bounding boxes for summary image
     24    // first, examine the PSF and SAT stars to set output image size:
     25    // - add stamp widths until we exceed output image width,
     26    // - then start a new row offset by max height
    2627    for (int i = 0; i < sources->n; i++) {
    2728
    … …  
    3132        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
    3233        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
    33         if (!keep) continue;
     34        if (!keep) {
     35            psTrace ("psphot", 4, "not plotting star %d: %x", i, source->mode);
     36            continue;
     37        }
    3438
    3539        // how does this subimage get placed into the output image?
    … …  
    5761    }
    5862
     63    if (NY == 0) {
     64        psWarning ("no PSF or SAT stars to plot? skipping.\n");
     65        return false;
     66    }
     67
    5968    // allocate output image
    6069    psImage *outpos = psImageAlloc (NX, NY, PS_TYPE_F32);
    6170    psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32);
     71
     72    psImageInit (outpos, 0.0);
     73    psImageInit (outsub, 0.0);
    6274
    6375    int Xo = 0;                         // starting corner of next box
    … …  
    142154    psFree (outpos);
    143155    psFree (outsub);
    144     return (sources);
     156    return true;
    145157}
    146158
Note: See TracChangeset for help on using the changeset viewer.