IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6481


Ignore:
Timestamp:
Feb 23, 2006, 6:16:11 PM (20 years ago)
Author:
eugene
Message:

fixed some errors with double sources, added second-pass linear PSF fit, multiple-depths

Location:
trunk/psphot/src
Files:
1 added
21 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/Makefile.am

    r6427 r6481  
    3737        psphotRadiusChecks.c    \
    3838        psphotImageMedian.c     \
     39        psphotSkyReplace.c      \
    3940        psLine.c                \
    4041        psSparse.c              \

  • trunk/psphot/src/pmModelFitSet.c

    r6379 r6481  
    4242    model = params->data.F32[0];
    4343
    44     PAR[0] = 0;
     44    PAR[0] = params->data.F32[0];
    4545    for (int i = 0; i < nSrc; i++) {
    4646        for (int n = 1; n < nPar; n++) {
    … …  
    5353        }
    5454    }
     55    deriv->data.F32[0] = dPAR[0]*2.0;
    5556    return (model);
    5657}

  • trunk/psphot/src/pmSourceFitSet.c

    r6441 r6481  
    7474    int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
    7575
     76    // define parameter vectors for source set
    7677    psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    7778    psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    78     psVector *paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL;
     79
     80    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
     81
     82    // define limits for single-source model
     83    psVector *oneDelta;
     84    psVector *oneParMin;
     85    psVector *oneParMax;
     86    modelLimits (&oneDelta, &oneParMin, &oneParMax);
     87
     88    psMinConstrain *constrain = psMinConstrainAlloc();
     89    constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
     90    constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
     91    constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
     92    constrain->paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL;
    7993
    8094    // all but the sky are allowed to vary independently (subject to PSF)
     95    // set the parameters from the multiple models
     96    // also, set limits based on single-source limits
    8197    params->data.F32[0] = model->params->data.F32[0];
     98    constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0];
     99    constrain->paramMin->data.F32[0]   = oneParMin->data.F32[0];
     100    constrain->paramMax->data.F32[0]   = oneParMax->data.F32[0];
     101    constrain->paramMask->data.U8[0]   = 0;
    82102    for (int i = 0; i < nSrc; i++) {
    83103      model = modelSet->data[i];
    … …  
    85105        params->data.F32[i*nPar + n] = model->params->data.F32[n];
    86106        dparams->data.F32[i*nPar + n] = model->dparams->data.F32[n];
     107        constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n];
     108        constrain->paramMin->data.F32[i*nPar + n]   = oneParMin->data.F32[n];
     109        constrain->paramMax->data.F32[i*nPar + n]   = oneParMax->data.F32[n];
    87110        if (PSF) {
    88             paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1;
     111            constrain->paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1;
    89112        }
    90113      }
     114    }
     115    psFree (oneDelta);
     116    psFree (oneParMin);
     117    psFree (oneParMax);
     118
     119    if (psTraceGetLevel(__func__) >= 5) {
     120        for (int i = 0; i < params->n; i++) {
     121            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]);
     122        }
    91123    }
    92124
    … …  
    126158    }
    127159    if (nPix <  nParams + 1) {
    128         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    129         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     160        psTrace (__func__, 4, "insufficient valid pixels\n");
     161        psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__);
    130162        model->status = PM_MODEL_BADARGS;
    131163        psFree (x);
    … …  
    134166        psFree (params);
    135167        psFree (dparams);
    136         psFree (paramMask);
     168        psFree(constrain->paramMask);
     169        psFree(constrain->paramMin);
     170        psFree(constrain->paramMax);
     171        psFree(constrain->paramDelta);
     172        psFree(constrain);
    137173        return(false);
    138174    }
    … …  
    143179    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    144180                            PM_SOURCE_FIT_MODEL_TOLERANCE);
    145     psMinConstrain *constrain = psMinConstrainAlloc();
    146     constrain->paramMask = paramMask;
    147 
    148     // Set the parameter range checks
    149     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    150     modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
    151181
    152182    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    153183
    154     psTrace (".pmObjects.pmSourceFitSet", 5, "fitting function\n");
     184    psTrace (__func__, 5, "fitting function\n");
    155185    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet);
    156186
    157187    // parameter errors from the covariance matrix
    158188    for (int i = 0; i < dparams->n; i++) {
    159         if ((paramMask != NULL) && paramMask->data.U8[i])
     189        if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i])
    160190            continue;
    161191        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    … …  
    163193
    164194    // get the Gauss-Newton distance for fixed model parameters
    165     if (paramMask != NULL) {
     195    if (constrain->paramMask != NULL) {
    166196        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    167197        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet);
    168198        for (int i = 0; i < dparams->n; i++) {
    169             if (!paramMask->data.U8[i])
     199            if (!constrain->paramMask->data.U8[i])
    170200                continue;
    171201            dparams->data.F32[i] = delta->data.F64[i];
    … …  
    220250
    221251    rc = (onPic && fitStatus);
    222     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     252    psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);
     253    psTrace(__func__, 3, "---- %s end : status %d ----\n", __func__, rc);
    223254    return(rc);
    224255}

  • trunk/psphot/src/psModulesUtils.c

    r6379 r6481  
    144144    int xIs, xJs, yIs, yJs;
    145145    int xIe, yIe;
    146     float flux;
     146    float flux, wt;
    147147
    148148    psImage *Pi = Mi->pixels;
    149149    psImage *Pj = Mj->pixels;
     150
     151    psImage *Wi = Mi->weight;
    150152
    151153    psImage *Ti = Mi->mask;
    … …  
    166168    yIe = Ye - Pi->row0;
    167169
     170    // note that this is addressing the same image pixels,
     171    // though only if both are source not model images
    168172    flux = 0;
    169173    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    … …  
    171175            if (Ti->data.U8[yi][xi]) continue;
    172176            if (Tj->data.U8[yj][xj]) continue;
    173             flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj];
     177            wt = Wi->data.F32[yi][xi];
     178            if (wt > 0) {
     179                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
     180            }
     181        }
     182    }
     183    return (flux);
     184}
     185
     186float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) {
     187
     188    int Xs, Xe, Ys, Ye;
     189    int xi, xj, yi, yj;
     190    int xIs, xJs, yIs, yJs;
     191    int xIe, yIe;
     192    float flux, wt;
     193
     194    psImage *Pi = Mi->pixels;
     195    psImage *Pj = Mj->pixels;
     196
     197    psImage *Wi = Mi->weight;
     198
     199    psImage *Ti = Mi->mask;
     200    psImage *Tj = Mj->mask;
     201
     202    Xs = PS_MAX (Pi->col0, Pj->col0);
     203    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
     204   
     205    Ys = PS_MAX (Pi->row0, Pj->row0);
     206    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
     207   
     208    xIs = Xs - Pi->col0;
     209    xJs = Xs - Pj->col0;
     210    yIs = Ys - Pi->row0;
     211    yJs = Ys - Pj->row0;
     212
     213    xIe = Xe - Pi->col0;
     214    yIe = Ye - Pi->row0;
     215
     216    // note that this is addressing the same image pixels,
     217    // though only if both are source not model images
     218    flux = 0;
     219    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
     220        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
     221            if (Ti->data.U8[yi][xi]) continue;
     222            if (Tj->data.U8[yj][xj]) continue;
     223            wt = Wi->data.F32[yi][xi];
     224            if (wt > 0) {
     225                flux += 1.0 / wt;
     226            }
    174227        }
    175228    }
    … …  
    250303    return NULL;
    251304}
     305
     306bool psRegionIsNaN (psRegion region) {
     307
     308    if (!isfinite(region.x0)) return true;
     309    if (!isfinite(region.x1)) return true;
     310    if (!isfinite(region.y0)) return true;
     311    if (!isfinite(region.y1)) return true;
     312    return false;
     313}

  • trunk/psphot/src/psModulesUtils.h

    r6311 r6481  
    4444// unify with paul's image/header/metadata functions
    4545psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
    46 pmModel *pmModelSelect (pmSource *source);
     46pmModel     *pmModelSelect (pmSource *source);
    4747
    48 ppConfig       *ppConfigAlloc (void);
    49 ppFile         *ppFileAlloc (void);
    50 psMetadata    *pmReadoutGetHeader (pmReadout *readout);
     48ppConfig    *ppConfigAlloc (void);
     49ppFile      *ppFileAlloc (void);
     50psMetadata  *pmReadoutGetHeader (pmReadout *readout);
     51
     52bool psRegionIsNaN (psRegion region);
    5153
    5254# endif

  • trunk/psphot/src/psSparse.c

    r6379 r6481  
    3434    sparse->Bfj->data.F32[2] = B->data.F32[2];
    3535
    36     x = psSparseSolve (x, sparse, 0);
    37     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
    38 
    39     x = psSparseSolve (x, sparse, 1);
    40     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
    41 
    42     x = psSparseSolve (x, sparse, 2);
    43     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
    44 
    45     x = psSparseSolve (x, sparse, 3);
    46     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
    47 
    48     x = psSparseSolve (x, sparse, 4);
     36    psSparseConstraint constraint;
     37    constraint.paramMin   = -1e8;
     38    constraint.paramMax   = +1e8;
     39    constraint.paramDelta = +1e8;
     40
     41    x = psSparseSolve (x, constraint, sparse, 0);
     42    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     43
     44    x = psSparseSolve (x, constraint, sparse, 1);
     45    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     46
     47    x = psSparseSolve (x, constraint, sparse, 2);
     48    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     49
     50    x = psSparseSolve (x, constraint, sparse, 3);
     51    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     52
     53    x = psSparseSolve (x, constraint, sparse, 4);
    4954    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
    5055    return;
    … …  
    166171}
    167172
    168 psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
     173psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter) {
    169174
    170175    psF32 dG;
    … …  
    182187        for (int i = 0; i < dQ->n; i++) {
    183188            dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
     189            if (fabs (dG) > constraint.paramDelta) {
     190                if (dG > 0) {
     191                    dG = +constraint.paramDelta;
     192                } else {
     193                    dG = -constraint.paramDelta;
     194                }
     195            }
    184196            guess->data.F32[i] -= dG;
     197            guess->data.F32[i] = PS_MAX (guess->data.F32[i], constraint.paramMin);
     198            guess->data.F32[i] = PS_MIN (guess->data.F32[i], constraint.paramMax);
    185199        }
    186200    }
    … …  
    220234}
    221235
    222 

  • trunk/psphot/src/psSparse.h

    r5654 r6481  
     1
     2typedef struct {
     3    double paramDelta;
     4    double paramMin;
     5    double paramMax;
     6} psSparseConstraint;
    17
    28typedef struct {
    … …  
    1521void psSparseVectorElement (psSparse *sparse, int i, float value);
    1622psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector);
    17 psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter);
     23psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter);
    1824psSparse *psSparseAlloc (int Nrows, int Nelem);

  • trunk/psphot/src/psphot.c

    r6379 r6481  
    2626    psFree (config);
    2727    psphotCleanup ();
    28 
     28 
    2929    exit (0);
    3030}

  • trunk/psphot/src/psphot.h

    r6427 r6481  
    4242
    4343// optional object analysis steps
    44 bool            psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf);
     44bool            psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final);
    4545bool            psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf);
    4646bool            psphotReplaceUnfit (psArray *sources);
    … …  
    5656pmModel        *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);
    5757float           pmSourceCrossProduct (pmSource *Mi, pmSource *Mj);
     58float           pmSourceCrossWeight (pmSource *Mi, pmSource *Mj);
    5859psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
    5960void            psphotModelGroupInit (void);
    … …  
    6364void            psphotTestArguments (int *argc, char **argv);
    6465bool            pmCellSetMask (pmCell *cell, psMetadata *recipe);
     66bool            psphotBackgroundNames (psMetadata *arguments);
     67bool            psphotSkyReplace (pmReadout *readout, psImage *background);
    6568
    6669// functions to set the correct source pixels
    6770bool            psphotInitRadiusPSF (psMetadata *config, pmModelType type);
    6871bool            psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model);
     72bool            psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR);
    6973bool            psphotInitRadiusEXT (psMetadata *config, pmModelType type);
    7074bool            psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model);

  • trunk/psphot/src/psphotArguments.c

    r6379 r6481  
    4242        psArgumentRemove (N, argc, argv);
    4343        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "RESID_IMAGE", PS_META_REPLACE, "", argv[N]);
     44        psArgumentRemove (N, argc, argv);
     45    }
     46
     47    // background image (output) is used by psphotOutput
     48    if ((N = psArgumentGet (*argc, argv, "-background"))) {
     49        psArgumentRemove (N, argc, argv);
     50        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_IMAGE", PS_META_REPLACE, "", argv[N]);
     51        psArgumentRemove (N, argc, argv);
     52    }
     53
     54    // background model (output) is used by psphotOutput
     55    if ((N = psArgumentGet (*argc, argv, "-backmodel"))) {
     56        psArgumentRemove (N, argc, argv);
     57        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_MODEL", PS_META_REPLACE, "", argv[N]);
    4458        psArgumentRemove (N, argc, argv);
    4559    }
    … …  
    89103        psArgumentRemove (N, argc, argv);
    90104        psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_INPUT_FILE", PS_META_REPLACE, "", argv[N]);
     105        psArgumentRemove (N, argc, argv);
     106    }
     107
     108    // PSF : optional psf file to be loaded
     109    if ((N = psArgumentGet (*argc, argv, "-psfout"))) {
     110        psArgumentRemove (N, argc, argv);
     111        psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_OUTPUT_FILE", PS_META_REPLACE, "", argv[N]);
    91112        psArgumentRemove (N, argc, argv);
    92113    }

  • trunk/psphot/src/psphotBasicDeblend.c

    r6427 r6481  
    8686        threshold = FRACTION * source->moments->Peak;
    8787        // threshold is no less than NSIGMA above the local median sigma?
    88         threshold = PS_MAX (threshold, NSIGMA*source->moments->dSky);
     88        threshold = PS_MAX (threshold, NSIGMA*sqrt(source->moments->dSky));
    8989
    9090        // generate a basic contour (set of x,y coordinates at-or-below flux level)

  • trunk/psphot/src/psphotBlendFit.c

    r6427 r6481  
    33// XXX I don't like this name
    44bool psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) {
     5
     6    bool status;
    57
    68    psTimerStart ("psphot");
    … …  
    911    sources = psArraySort (sources, psphotSortBySN);
    1012   
     13    // S/N limit to perform full non-linear fits
     14    float FIT_SN_LIM = psMetadataLookupF32 (&status, config, "FULL_FIT_SN_LIM");
     15
    1116    psphotInitLimitsPSF (config);
    1217    psphotInitLimitsEXT (config);
    1318    psphotInitRadiusPSF (config, psf->type);
    1419
     20    // option to limit analysis to a specific region
     21    char *region = psMetadataLookupStr (&status, config, "ANALYSIS_REGION");
     22    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     23    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
     24
    1525    for (int i = 0; i < sources->n; i++) {
    16 
    1726        pmSource *source = sources->data[i];
    1827
    … …  
    2231        if (source->type == PM_SOURCE_SATURATED) continue;
    2332
     33        // limit selection to some SN limit
     34        if (source->moments == NULL) continue;
     35        if (source->moments->SN < FIT_SN_LIM) continue;
     36
    2437        // if model is NULL, we don't have a starting guess
    2538        if (source->modelPSF == NULL) continue;
     39
     40        if (source->moments->x < AnalysisRegion.x0) continue;
     41        if (source->moments->y < AnalysisRegion.y0) continue;
     42        if (source->moments->x > AnalysisRegion.x1) continue;
     43        if (source->moments->y > AnalysisRegion.y1) continue;
    2644
    2745        // replace object in image
    … …  
    2947
    3048        psTrace ("psphot.blend", 5, "trying source at %f, %f\n", source->moments->x, source->moments->y);
    31 
     49       
    3250        // try fitting PSFs, then try extended sources
    3351        if (psphotFitBlend (readout, source)) continue;

  • trunk/psphot/src/psphotChoosePSF.c

    r6427 r6481  
    1111
    1212    psTimerStart ("psphot");
     13
     14    // check if a PSF model is supplied by the user
     15    char *psfFile = psMetadataLookupStr (&status, config, "PSF_INPUT_FILE");
     16    if (status) {
     17        int Nfail;
     18        psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, psfFile, FALSE);
     19        pmPSF *psf = pmPSFfromMD (psfData);
     20        return psf;
     21    }
    1322
    1423    // examine PSF sources in S/N order (brightest first)

  • trunk/psphot/src/psphotEnsemblePSF.c

    r6427 r6481  
    11# include "psphot.h"
     2bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);
    23
    34// 2006.02.07 : no leaks!
    45// fit all reasonable sources with the linear PSF model
    5 bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) {
     6bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final) {
    67
    78    bool  status;
    … …  
    910    float y;
    1011    float f;
     12    float r;
     13
     14    // psRegion allArray = psRegionSet (0, 0, 0, 0);
    1115
    1216    psTimerStart ("psphot");
    … …  
    2933
    3034    // option to limit analysis to a specific region
    31     bool UseAnalysisRegion = false;
    32     psRegion AnalysisRegion = {0, 0, 0, 0};
    3335    char *region = psMetadataLookupStr (&status, config, "ANALYSIS_REGION");
    34     if (status) {
    35         UseAnalysisRegion = true;
    36         AnalysisRegion = psRegionFromString (region);
    37         psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n",
    38                   AnalysisRegion.x0, AnalysisRegion.y0,
    39                   AnalysisRegion.x1, AnalysisRegion.y1);
    40     }
     36    psRegion AnalysisRegion = psRegionFromString (region);
     37    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
     38    // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     39    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
    4140
    4241    for (int i = 0; i < sources->n; i++) {
    … …  
    4544        // skip non-astronomical objects (very likely defects)
    4645        // XXX EAM : should we try these anyway?
    47         if (inSource->mode &  PM_SOURCE_BLEND) continue;
    4846        if (inSource->type == PM_SOURCE_DEFECT) continue;
    4947        if (inSource->type == PM_SOURCE_SATURATED) continue;
    50 
    51         if (UseAnalysisRegion) {
    52             if (inSource->moments->x < AnalysisRegion.x0) continue;
    53             if (inSource->moments->y < AnalysisRegion.y0) continue;
    54             if (inSource->moments->x > AnalysisRegion.x1) continue;
    55             if (inSource->moments->y > AnalysisRegion.y1) continue;
    56         }
     48        if (final) {
     49            if (inSource->mode &  PM_SOURCE_SUBTRACTED) continue;
     50        } else {
     51            if (inSource->mode &  PM_SOURCE_BLEND) continue;
     52        }
     53
     54        if (inSource->moments->x < AnalysisRegion.x0) continue;
     55        if (inSource->moments->y < AnalysisRegion.y0) continue;
     56        if (inSource->moments->x > AnalysisRegion.x1) continue;
     57        if (inSource->moments->y > AnalysisRegion.y1) continue;
    5758
    5859        pmSource *otSource = pmSourceAlloc ();
    … …  
    7374            // peak-up on peak (for non-sat objects)
    7475
    75             int ix = inSource->peak->x;
    76             int iy = inSource->peak->y;
     76            // ix,iy must land on inSource->pixels
     77            int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
     78            int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));
    7779
    7880            psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
    … …  
    104106        otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
    105107        otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
     108        otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);
    106109
    107110        // build the model image
    … …  
    117120        psImageKeepCircle (mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
    118121        pmSourceAddModel (flux, mask, model, false, false);
     122
     123        // calculate nDOF (nPix - 1)
     124        // int Nmaskpix = psImageCountPixelMask (mask, allArray, PSPHOT_MASK_SATURATED);
     125        // model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;
    119126
    120127        // save source in list
    … …  
    128135    // fill out the sparse matrix
    129136    psSparse *sparse = psSparseAlloc (models->n, 100);
     137    psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
    130138    for (int i = 0; i < models->n; i++) {
    131139        int N = index->data.U32[i];
    … …  
    133141        pmSource *Mi = models->data[i];
    134142
    135         // diagonal element (auto-cross-product)
    136         f = pmSourceCrossProduct (Mi, Mi);
    137         psSparseMatrixElement (sparse, i, i, f);
     143        // scale by diagonal element (auto-cross-product)
     144        r = pmSourceCrossProduct (Mi, Mi);
     145        weight->data.F32[i] = r;
     146
     147        psSparseMatrixElement (sparse, i, i, 1.0);
    138148
    139149        // find the image x model value
    140150        f = pmSourceCrossProduct (Fi, Mi);
    141         psSparseVectorElement (sparse, i, f);
     151        psSparseVectorElement (sparse, i, f / r);
    142152
    143153        // loop over all other stars following this one
    … …  
    153163            // got an overlap; calculate cross-product and add to output array
    154164            f = pmSourceCrossProduct (Mi, Mj);
    155             psSparseMatrixElement (sparse, j, i, f);
     165            psSparseMatrixElement (sparse, j, i, f / r);
    156166        }
    157167    }
    158168    psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
     169
     170    psSparseConstraint constraint;
     171    constraint.paramMin   = 0;
     172    constraint.paramMax   = 1e8;
     173    constraint.paramDelta = 1e8;
    159174
    160175    // solve for normalization terms (need include local sky?)
    161176    psSparseResort (sparse);
    162     psVector *norm = psSparseSolve (NULL, sparse, 3);
     177    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
    163178
    164179    // adjust models, set sources and subtract
    … …  
    168183        pmSource *Mi = models->data[i];
    169184
     185        // if we already have a PSF model, free it.
     186        psFree (Fi->modelPSF);
     187
    170188        // need to increment counter so we can free models here and sources above
    171189        Fi->modelPSF = psMemCopy (Mi->modelPSF);
    … …  
    173191        // assign linearly-fitted normalization
    174192        Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
     193        Fi->modelPSF->dparams->data.F32[1] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
     194        // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
    175195
    176196        // subtract object
    177197        pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
    178198        Fi->mode |= PM_SOURCE_SUBTRACTED;
    179         Fi->mode |= PM_SOURCE_TEMPSUB;
    180     }
     199        if (!final) Fi->mode |= PM_SOURCE_TEMPSUB;
     200    }
     201
     202    // measure chisq for each source
     203    for (int i = 0; final && (i < models->n); i++) {
     204        int N = index->data.U32[i];
     205        pmSource *Fi = sources->data[N];
     206        pmModel *model = Fi->modelPSF;
     207
     208        x = model->params->data.F32[2];
     209        y = model->params->data.F32[3];
     210
     211        psImageKeepCircle (Fi->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
     212        pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
     213        psImageKeepCircle (Fi->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
     214    }
     215
    181216    psFree (index);
    182217    psFree (sparse);
    183218    psFree (models);
    184219    psFree (norm);
     220    psFree (weight);
    185221
    186222    psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
    187223    return true;
    188224}
     225
     226bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight) {
     227
     228    double dC = 0.0;
     229    int Npix = 0;
     230    for (int j = 0; j < image->numRows; j++) {
     231        for (int i = 0; i < image->numCols; i++) {
     232            if (mask->data.U8[j][i]) continue;
     233            if (weight->data.F32[j][i] <= 0) continue;
     234            dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
     235            Npix ++;
     236        }
     237    }
     238    model->nDOF = Npix - 1;
     239    model->chisq = dC;
     240
     241    return (true);
     242}

  • trunk/psphot/src/psphotFindPeaks.c

    r6427 r6481  
    1010
    1111    // smooth the image and weight map
    12 
    13     psphotSaveImage (NULL, readout->image, "image.fits");
    14 
    1512    psTimerStart ("psphot");
    1613
    … …  
    2219    psLogMsg ("psphot", 4, "smooth image: %f sec\n", psTimerMark ("psphot"));
    2320
    24     psphotSaveImage (NULL, smooth_im, "smooth_im.fits");
    25 
    2621    psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
    2722    psImageSmooth (smooth_wt, SIGMA, NSIGMA);
    2823    psLogMsg ("psphot", 4, "smooth weight: %f sec\n", psTimerMark ("psphot"));
    29 
    30     psphotSaveImage (NULL, smooth_wt, "smooth_wt.fits");
    3124
    3225    psTimerStart ("psphot");

  • trunk/psphot/src/psphotImageMedian.c

    r6441 r6481  
    88// use no more than MAX_SAMPLE_PIXELS pixels for each median box
    99static int MAX_SAMPLE_PIXELS;
     10
     11static char *backImage = NULL;
     12static char *backModel = NULL;
     13
     14bool psphotBackgroundNames (psMetadata *arguments) {
     15
     16    bool status;
     17
     18    backImage = psMetadataLookupStr (&status, arguments, "BACKGROUND_IMAGE");
     19    backModel = psMetadataLookupStr (&status, arguments, "BACKGROUND_MODEL");
     20
     21    return true;
     22}
    1023
    1124// generate the median in NxN boxes, clipping heavily
    … …  
    211224        }
    212225    }
     226
     227    // optionally save background
     228    if (backImage != NULL) psphotSaveImage (NULL, background, backImage);
     229    if (backModel != NULL) psphotSaveImage (NULL, model, backModel);
     230
    213231    psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
    214232    psFree (rnd);

  • trunk/psphot/src/psphotModelTest.c

    r6379 r6481  
    2626    // in psfMode, psf sets the model type
    2727    if (psfMode) {
    28         char *psfFile = psMetadataLookupStr (&status, arguments, "PSF_INPUT_FILE");
     28        char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE");
    2929        if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
    3030        psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, psfFile, FALSE);
    … …  
    157157    if (status) {
    158158        status = psphotFitSet (source, model, fitset, psfMode);
    159         return status;
     159        exit (0);
    160160    }
    161161

  • trunk/psphot/src/psphotOutput.c

    r6379 r6481  
    1111static char *extNameKey = NULL;
    1212static char *extRoot    = NULL;
     13static char *psfFile    = NULL;
     14static char *psfSample  = NULL;
    1315
    1416void psphotOutputCleanup () {
    … …  
    2123    psFree (extNameKey);
    2224    psFree (extRoot);
     25    psFree (psfFile);
     26    psFree (psfSample);
    2327    return;
    2428}
    … …  
    3438    outputMode      = psMetadataLookupStr (&status, config->recipe, "OUTPUT_MODE");
    3539    outputFormat    = psMetadataLookupStr (&status, config->recipe, "OUTPUT_FORMAT");
     40
     41    psfFile         = psMetadataLookupStr (&status, config->recipe, "PSF_OUTPUT_FILE");
     42    psfSample       = psMetadataLookupStr (&status, config->recipe, "PSF_SAMPLE_FILE");
    3643
    3744    // rules to construct output names
    … …  
    7582    }
    7683
     84    // register background image-file names
     85    psphotBackgroundNames (config->arguments);
     86
    7787    // save so freeing config will not drop these references
    7888    psMemCopy (outputRoot);
    … …  
    8393    psMemCopy (extNameKey);
    8494    psMemCopy (extRoot);
     95    psMemCopy (psfFile);
     96    psMemCopy (psfSample);
    8597
    8698    return;
    … …  
    165177    char *outputFile;
    166178
    167     psTimerStart ("psphot");
    168 
    169179    psMetadata *header = pmReadoutGetHeader (readout);
    170180
    … …  
    180190    }
    181191
    182     char *psfFile    = psMetadataLookupStr (&status, arguments, "PSF_OUTPUT_FILE");
    183     char *psfSample  = psMetadataLookupStr (&status, arguments, "PSF_SAMPLE_FILE");
    184192    char *residImage = psMetadataLookupStr (&status, arguments, "RESID_IMAGE");
    185193
    … …  
    236244    psF32 *PAR, *dPAR;
    237245    float dmag, apResid;
     246
     247    psTimerStart ("string");
    238248
    239249    psLine *line = psLineAlloc (104);  // 104 is dophot-defined line length
    … …  
    274284    }
    275285    fclose (f);
     286    psFree (line);
     287    fprintf (stderr, "%f seconds for %d objects with psLine\n", psTimerMark ("string"), (int)sources->n);
     288
     289    psTimerStart ("string");
     290
     291    f = fopen ("test.obj", "w");
     292    if (f == NULL) {
     293        psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", "test.obj");
     294        return false;
     295    }
     296
     297    char *string;
     298    // write sources with models
     299    for (int i = 0; i < sources->n; i++) {
     300        pmSource *source = (pmSource *) sources->data[i];
     301        pmModel *model = pmModelSelect (source);
     302        if (model == NULL) continue;
     303
     304        PAR = model->params->data.F32;
     305        dPAR = model->dparams->data.F32;
     306
     307        dmag = dPAR[1] / PAR[1];
     308        type = pmSourceDophotType (source);
     309        apResid = source->apMag - source->fitMag;
     310
     311        string = NULL;
     312        psStringAppend (&string, "%3d",   type);
     313        psStringAppend (&string, "%8.2f", PAR[2]);
     314        psStringAppend (&string, "%8.2f", PAR[3]);
     315        psStringAppend (&string, "%8.3f", source->fitMag);
     316        psStringAppend (&string, "%6.3f", dmag);
     317        psStringAppend (&string, "%9.2f", PAR[0]);
     318        psStringAppend (&string, "%9.3f", PAR[4]);
     319        psStringAppend (&string, "%9.3f", PAR[5]);
     320        psStringAppend (&string, "%7.2f", PAR[6]);
     321        psStringAppend (&string, "%8.3f", 99.999);
     322        psStringAppend (&string, "%8.3f", source->apMag);
     323        psStringAppend (&string, "%8.2f\n", apResid);
     324        fwrite (string, 1, strlen(string), f);
     325        psFree (string);
     326    }
     327    fclose (f);
     328    fprintf (stderr, "%f seconds for %d objects with psString\n", psTimerMark ("string"), (int)sources->n);
     329
    276330    return true;
    277331}
    … …  
    501555    }
    502556
    503    // write sources with models first
     557    // write sources with models first
    504558    for (i = 0; i < sources->n; i++) {
    505559        pmSource *source = (pmSource *) sources->data[i];
    … …  
    525579            fprintf (f, "%9.6f ", dPAR[j]);
    526580        }
    527         fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
     581        fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
    528582                 source[0].apMag, source[0].type, source[0].mode,
    529583                 log10(model[0].chisq/model[0].nDOF),
    … …  
    617671           
    618672        if (source->moments == NULL) {
    619           moment = empty;
     673            moment = empty;
    620674        } else {
    621           moment = source->moments;
     675            moment = source->moments;
    622676        }
    623677

  • trunk/psphot/src/psphotRadiusChecks.c

    r6427 r6481  
    2424    // set the fit radius based on the object flux limit and the model
    2525    model->radius = modelRadiusPSF (model->params, PSF_FIT_NSIGMA*moments->dSky) + PSF_FIT_PADDING;
     26    if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius");
     27       
     28    if (source->mode &  PM_SOURCE_SATSTAR) {
     29        model->radius *= 2;
     30    }
     31
     32    bool status = psphotRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radius);
     33    return status;
     34}
     35
     36bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR) {
     37
     38    pmMoments *moments = source->moments;
     39    if (moments == NULL) return false;
     40
     41    // set the fit radius based on the object flux limit and the model
     42    model->radius = modelRadiusPSF (model->params, PSF_FIT_NSIGMA*moments->dSky) + dR + PSF_FIT_PADDING;
    2643    if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius");
    2744       

  • trunk/psphot/src/psphotReadout.c

    r6442 r6481  
    1717
    1818    // construct sources and measure basic stats
     19    // limit moments analysis by S/N?
    1920    sources = psphotSourceStats (readout, config, peaks);
    2021
    2122    // classify sources based on moments, brightness
     23    // faint sources not classified?
    2224    psphotRoughClass (sources, config);
    2325
    … …  
    2830    psf = psphotChoosePSF (config, sources);
    2931
    30     // select analysis method
    31     char *FITMODE = psMetadataLookupStr (&status, config, "FITMODE");
    32     if (!strcasecmp(FITMODE, "ENSEMBLE")) {
    33         psphotEnsemblePSF (readout, config, sources, psf);
    34     }
     32    // linear PSF fit to peaks
     33    psphotEnsemblePSF (readout, config, sources, psf, FALSE);
    3534
    36     if (!strcasecmp(FITMODE, "BLEND")) {
    37         psphotEnsemblePSF (readout, config, sources, psf);
    38         psphotBlendFit (readout, config, sources, psf);
    39         psphotReplaceUnfit (sources);
    40         psphotApResid (readout, sources, config, psf);
    41     }
     35    // non-linear PSF and EXT fit to brighter sources
     36    psphotBlendFit (readout, config, sources, psf);
    4237
     38    // replace fitted sources
     39    psphotReplaceUnfit (sources);
     40
     41    // find remaining peaks after first source subtraction pass
     42    // psphotFindPeaks ();
     43
     44    // linear PSF fit to remaining peaks
     45    psphotEnsemblePSF (readout, config, sources, psf, TRUE);
     46
     47    // measure aperture photometry corrections
     48    psphotApResid (readout, sources, config, psf);
     49
     50    // calculate source magnitudes
    4351    psphotMagnitudes (config, sources, psf);
    4452
     53    // update the header with stats results
    4554    psMetadata *header = pmReadoutGetHeader (readout);
    4655    psphotUpdateHeader (header, config);
    4756
    48     // psphotSkyReplace (readout, skymodel);
     57    // XXX make this an option?
     58    // replace background in residual image
     59    psphotSkyReplace (readout, skymodel);
    4960
    5061    // need to do something with the sources, psf, and sky
    … …  
    6677// XXX Deprecate or allow these versions?
    6778# if (0)
     79    // select analysis method
     80    char *FITMODE = psMetadataLookupStr (&status, config, "FITMODE");
     81    if (!strcasecmp(FITMODE, "ENSEMBLE")) {
     82        psphotEnsemblePSF (readout, config, sources, psf);
     83    }
     84
    6885    if (!strcasecmp(FITMODE, "FULL")) {
    6986        psphotEnsemblePSF (readout, config, sources, psf);

  • trunk/psphot/src/psphotSourceFits.c

    r6427 r6481  
    126126    psTrace ("psphot.blend", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[2], ONE->params->data.F32[3]);
    127127
    128     source->modelPSF = (pmModel *) DBL->data[0];
    129     source->mode |=  PM_SOURCE_SUBTRACTED;
    130     source->mode &= ~PM_SOURCE_TEMPSUB;
    131 
     128    // drop old model, save new second model...
     129    psFree (source->modelPSF);
     130    source->modelPSF = psMemCopy (DBL->data[0]);
     131    source->mode    |= PM_SOURCE_SUBTRACTED;
     132    source->mode    |= PM_SOURCE_PAIR;
     133    source->mode    &= ~PM_SOURCE_TEMPSUB;
     134
     135    // copy most data from the primary source (modelEXT, blends stay NULL)
     136    pmSource *newSrc = pmSourceAlloc ();
     137    newSrc->peak     = psMemCopy (source->peak);
     138    newSrc->pixels   = psMemCopy (source->pixels);
     139    newSrc->weight   = psMemCopy (source->weight);
     140    newSrc->mask     = psMemCopy (source->mask);
     141    newSrc->moments  = psMemCopy (source->moments);
     142    newSrc->modelPSF = psMemCopy (DBL->data[1]);
     143    newSrc->type     = source->type;
     144    newSrc->mode     = source->mode;
     145    psArrayAdd (sources, 100, newSrc);
     146    psFree (newSrc);
     147    psFree (DBL);
    132148    return true;
    133149}
    … …  
    154170    // save the PSF model from the Ensemble fit
    155171    PSF = source->modelPSF;
    156     psphotCheckRadiusPSF (readout, source, PSF);
     172    psphotCheckRadiusPSFBlend (readout, source, PSF, 8.0);
    157173
    158174    modelSet = psArrayAlloc (2);
    … …  
    160176    DBL = pmModelCopy (PSF);
    161177    DBL->params->data.F32[1] *= 0.5;
    162     DBL->params->data.F32[2] += dx;
    163     DBL->params->data.F32[3] += dy;
     178    DBL->params->data.F32[2] = source->moments->x + dx;
     179    DBL->params->data.F32[3] = source->moments->y + dy;
    164180    modelSet->data[0] = DBL;
    165181
    166182    DBL = pmModelCopy (PSF);
    167183    DBL->params->data.F32[1] *= 0.5;
    168     DBL->params->data.F32[2] -= dx;
    169     DBL->params->data.F32[3] -= dy;
     184    DBL->params->data.F32[2] = source->moments->x - dx;
     185    DBL->params->data.F32[3] = source->moments->y - dy;
    170186    modelSet->data[1] = DBL;
    171187
    172     x = PSF->params->data.F32[2];
    173     y = PSF->params->data.F32[3];
     188    x = source->moments->x;
     189    y = source->moments->y;
    174190
    175191    // fit EXT (not PSF) model (set/unset the pixel mask)
    … …  
    214230    // save the PSF model from the Ensemble fit
    215231    pmModel *PSF = pmModelCopy (source->modelPSF);
    216 
    217     // extend source radius as needed
    218     psphotCheckRadiusPSF (readout, source, PSF);
    219232
    220233    x = PSF->params->data.F32[2];
    … …  
    229242    psArrayAdd (sourceSet, 16, source);
    230243
    231     // counter to track the blend elements used in the fit
    232     psVector *index  = psVectorAlloc (source->blends->n + 1, PS_TYPE_S32);
    233     index->data.S32[0] = -1; // first element corresponds to the primary source
    234 
     244    // we need to include all blends in the fit (unless primary is saturated?)
     245    dR = 0;
    235246    for (int i = 0; i < source->blends->n; i++) {
    236247        pmSource *blend = source->blends->data[i];
    237248
    238         // is this object close enough to include in the fit
    239         // XXX this test is bogus: need to think about this
    240         dR = hypot (blend->peak->x - x, blend->peak->y - y);
    241         if (dR > PSF->radius - 3) continue;
     249        // find the blend which is furthest from source
     250        dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
    242251
    243252        // create the model and guess parameters for this blend
    … …  
    256265
    257266        // add this blend to the list
    258         index->data.S32[modelSet->n] = i;
    259267        psArrayAdd (modelSet, 16, model);
    260268        psArrayAdd (sourceSet, 16, blend);
    261269    }
     270
     271    // extend source radius as needed
     272    psphotCheckRadiusPSFBlend (readout, source, PSF, dR);
    262273
    263274    // fit PSF model (set/unset the pixel mask)
    … …  
    273284        // if we skip this one, free the corresponding blend entry model
    274285        if (!psphotEvalPSF (src, model)) {
    275             int n = index->data.S32[i];
    276             pmSource *blend = source->blends->data[n];
     286            pmSource *blend = source->blends->data[i - 1];
    277287            psFree (blend->modelPSF);
    278288            blend->modelPSF = NULL;
    … …  
    285295        src->mode &= ~PM_SOURCE_TEMPSUB;
    286296    }
    287     psFree (index);
    288297    psFree (modelSet);
    289298    psFree (sourceSet);
Note: See TracChangeset for help on using the changeset viewer.