IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 16, 2006, 4:24:03 PM (20 years ago)
Author:
Paul Price
Message:

NOT_U8 and NOT_U16 macros were moved into psLib (and renamed with PS_ prefix). Updating files that use these macros so that they compile.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotSourceFits.c

    r9576 r9595  
    33// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
    44
    5 bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 
     5bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) {
    66
    77    float x, y, dR;
    … …  
    99    // if this source is not a possible blend, just fit as PSF
    1010    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
    11         bool status = psphotFitPSF (readout, source, psf);
    12         return status;
     11        bool status = psphotFitPSF (readout, source, psf);
     12        return status;
    1313    }
    1414    psTrace ("psphot", 5, "trying blend...\n");
    … …  
    2323    psArray *modelSet = psArrayAlloc (source->blends->n + 1);
    2424    psArrayAdd (modelSet, 16, PSF);
    25    
     25
    2626    psArray *sourceSet = psArrayAlloc (source->blends->n + 1);
    2727    psArrayAdd (sourceSet, 16, source);
    … …  
    3030    dR = 0;
    3131    for (int i = 0; i < source->blends->n; i++) {
    32         pmSource *blend = source->blends->data[i];
    33 
    34         // find the blend which is furthest from source
    35         dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
    36 
    37         // create the model and guess parameters for this blend
    38         pmModel *model = pmModelAlloc (PSF->type);
    39         for (int j = 0; j < model->params->n; j++) {
    40             model->params->data.F32[j] = PSF->params->data.F32[j];
    41             model->dparams->data.F32[j] = PSF->dparams->data.F32[j];
    42         }
    43 
    44         // XXX assume local sky is 0.0?
    45         model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
    46         if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit");
    47         model->params->data.F32[2] = blend->peak->x;
    48         model->params->data.F32[3] = blend->peak->y;
    49 
    50         // add this blend to the list
    51         psArrayAdd (modelSet, 16, model);
    52         psArrayAdd (sourceSet, 16, blend);
    53 
    54         // free to avoid double counting model
    55         psFree (model);
     32        pmSource *blend = source->blends->data[i];
     33
     34        // find the blend which is furthest from source
     35        dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
     36
     37        // create the model and guess parameters for this blend
     38        pmModel *model = pmModelAlloc (PSF->type);
     39        for (int j = 0; j < model->params->n; j++) {
     40            model->params->data.F32[j] = PSF->params->data.F32[j];
     41            model->dparams->data.F32[j] = PSF->dparams->data.F32[j];
     42        }
     43
     44        // XXX assume local sky is 0.0?
     45        model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
     46        if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit");
     47        model->params->data.F32[2] = blend->peak->x;
     48        model->params->data.F32[3] = blend->peak->y;
     49
     50        // add this blend to the list
     51        psArrayAdd (modelSet, 16, model);
     52        psArrayAdd (sourceSet, 16, blend);
     53
     54        // free to avoid double counting model
     55        psFree (model);
    5656    }
    5757
    … …  
    6262    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
    6363    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
    64     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
     64    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK));
    6565
    6666    // correct model chisq for flux trend
    … …  
    7070    // evaluate the blend objects, subtract if good, free otherwise
    7171    for (int i = 1; i < modelSet->n; i++) {
    72         pmSource *blend = sourceSet->data[i];
    73         pmModel *model  = modelSet->data[i];
    74 
    75         // correct model chisq for flux trend
    76         chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
    77         model->chisqNorm = model->chisq / chiTrend;
    78 
    79         // if this one failed, skip it
    80         if (!psphotEvalPSF (blend, model)) {
    81             psTrace ("psphot", 5, "failed on blend %d of %ld\n", i, modelSet->n);
    82             continue;
    83         }
    84 
    85         // otherwise, supply the resulting model to the corresponding blend
    86         psFree(blend->modelPSF);
    87         blend->modelPSF = psMemIncrRefCounter (model);
    88         psTrace ("psphot", 5, "fitted blend as PSF\n");
    89         pmModelSub (source->pixels, source->mask, model, false, false);
    90         blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    91         blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
     72        pmSource *blend = sourceSet->data[i];
     73        pmModel *model  = modelSet->data[i];
     74
     75        // correct model chisq for flux trend
     76        chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
     77        model->chisqNorm = model->chisq / chiTrend;
     78
     79        // if this one failed, skip it
     80        if (!psphotEvalPSF (blend, model)) {
     81            psTrace ("psphot", 5, "failed on blend %d of %ld\n", i, modelSet->n);
     82            continue;
     83        }
     84
     85        // otherwise, supply the resulting model to the corresponding blend
     86        psFree(blend->modelPSF);
     87        blend->modelPSF = psMemIncrRefCounter (model);
     88        psTrace ("psphot", 5, "fitted blend as PSF\n");
     89        pmModelSub (source->pixels, source->mask, model, false, false);
     90        blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
     91        blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
    9292    }
    9393    psFree (modelSet);
    … …  
    9696    // evaluate the primary object
    9797    if (!psphotEvalPSF (source, PSF)) {
    98         psTrace ("psphot", 5, "failed on blend 0 of %ld\n", modelSet->n);
    99         psFree (PSF);
    100         return false;
     98        psTrace ("psphot", 5, "failed on blend 0 of %ld\n", modelSet->n);
     99        psFree (PSF);
     100        return false;
    101101    }
    102102
    … …  
    110110}
    111111
    112 bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { 
     112bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) {
    113113
    114114    float x, y;
    … …  
    128128    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
    129129    pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF);
    130     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
    131    
     130    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK));
     131
    132132    // correct model chisq for flux trend
    133133    chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
    … …  
    136136    // does the PSF model succeed?
    137137    if (!psphotEvalPSF (source, PSF)) {
    138         psFree (PSF);
    139         return false;
     138        psFree (PSF);
     139        return false;
    140140    }
    141141
    … …  
    143143    pmModelSub (source->pixels, source->mask, PSF, false, false);
    144144
    145     // free old model, save new model 
     145    // free old model, save new model
    146146    psFree (source->modelPSF);
    147147    source->modelPSF = PSF;
    … …  
    173173}
    174174
    175 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { 
     175bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) {
    176176
    177177    bool okEXT, okDBL;
    … …  
    216216    ONE->chisqNorm = ONE->chisq / chiTrend;
    217217
    218     psFree (tmpSrc); 
     218    psFree (tmpSrc);
    219219
    220220    if (okEXT && okDBL) {
    221         psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
    222         // XXX EAM : a bogus bias: need to examine this better
    223         if (3*chiEXT > chiDBL) goto keepDBL;
    224         goto keepEXT;
     221        psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
     222        // XXX EAM : a bogus bias: need to examine this better
     223        if (3*chiEXT > chiDBL) goto keepDBL;
     224        goto keepEXT;
    225225    }
    226226
    … …  
    276276
    277277// fit a double PSF source to an extended blob
    278 psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { 
     278psArray *psphotFitDBL (pmReadout *readout, pmSource *source) {
    279279
    280280    float x, y, dx, dy;
    … …  
    284284    psEllipseMoments moments;
    285285    psArray *modelSet;
    286    
     286
    287287    // make a guess at the position of the two sources
    288288    moments.x2 = source->moments->Sx;
    … …  
    320320    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
    321321    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
    322     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
     322    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK));
    323323
    324324    return (modelSet);
    325325}
    326326
    327 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { 
     327pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) {
    328328
    329329    float x, y;
    330330
    331331    // use the source moments, etc to guess basic model parameters
    332     pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); 
     332    pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT);
    333333    // if (isnan(EXT->params->data.F32[1])) psAbort ("psphot", "nan in ext fit");
    334334
    … …  
    339339
    340340    if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) {
    341         psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy);
     341        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy);
    342342    }
    343343
    … …  
    345345    psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "OR", PM_MASK_MARK);
    346346    pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT);
    347     psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
     347    psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK));
    348348
    349349    return (EXT);
Note: See TracChangeset for help on using the changeset viewer.