IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13035


Ignore:
Timestamp:
Apr 25, 2007, 3:35:50 PM (19 years ago)
Author:
eugene
Message:

incorporating changes from eam_02_branch (pmSourceAdd/Sub, modified API for pmModelAdd, cached models)

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

Legend:

Unmodified
Added
Removed
  • trunk/psphot/src/Makefile.am

    r13007 r13035  
    3131        psphotChoosePSF.c       \
    3232        psphotGuessModels.c     \
    33         psphotEnsemblePSF.c     \
     33        psphotFitSourcesLinear.c \
    3434        psphotBlendFit.c        \
    3535        psphotReplaceUnfit.c    \
     
    4747        psphotModelTest.c       \
    4848        psphotFitSet.c          \
    49         psphotWeightBias.c      \
    5049        psphotSourceFreePixels.c \
    5150        psphotSummaryPlots.c     \
    52         psphotTestPSF.c          \
    5351        psphotMergeSources.c     \
    5452        psphotReadoutCleanup.c   \
     
    5856        psphotMosaicSubimage.c   \
    5957        psphotMakeResiduals.c    \
    60         psphotTestSourceOutput.c \
    6158        psphotAddNoise.c
    6259
  • trunk/psphot/src/psphot.h

    r13007 r13035  
    4141bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
    4242bool            psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
     43bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
    4344bool            psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
    4445bool            psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
  • trunk/psphot/src/psphotAddNoise.c

    r12922 r13035  
    6767        PAR[PM_PAR_SXY] = newshape.sxy;
    6868
    69         if (add) {
    70             pmModelAdd (source->weight, source->mask, model, false, false);
    71         } else {
    72             pmModelSub (source->weight, source->mask, model, false, false);
    73         }
     69        // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected
     70        pmSourceOp (source, PM_MODEL_OP_FULL, add);
    7471       
    7572        // restore original values
  • trunk/psphot/src/psphotBlendFit.c

    r13006 r13035  
    4343
    4444        // limit selection to some SN limit
     45        // XXX this should use peak?
    4546        if (source->moments == NULL) continue;
    4647        if (source->moments->SN < FIT_SN_LIM) continue;
    4748
     49        // XXX this should use peak?
    4850        if (source->moments->x < AnalysisRegion.x0) continue;
    4951        if (source->moments->y < AnalysisRegion.y0) continue;
     
    6567        // replace object in image
    6668        if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
    67             pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false);
     69            pmSourceAdd (source, PM_MODEL_OP_FULL);
    6870        }
    6971        Nfit ++;
    7072
    7173        // try fitting PSFs, then try extended sources
     74        // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
    7275        if (psphotFitBlend (readout, source, psf)) {
    7376            psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y);
     
    8487        Nfail ++;
    8588
    86         // re-subtract PSF for object, leave local sky
    87         pmModelSub (source->pixels, source->mask, source->modelPSF, false, false);
     89        // re-subtract the object, leave local sky
     90        pmSourceCacheModel (source);
     91        pmSourceSub (source, PM_MODEL_OP_FULL);
    8892        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
    8993        source->mode |= PM_SOURCE_MODE_TEMPSUB;
  • trunk/psphot/src/psphotChoosePSF.c

    r13011 r13035  
    239239
    240240            // set the mask and subtract the PSF model
    241             psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
    242             pmModelSub (source->pixels, source->mask, source->modelPSF, false, false);
    243             psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
     241            // XXX should we be using maskObj? should we be unsetting the mask?
     242            // use pmModelSub because modelFlux has not been generated
     243            assert (source->maskObj);
     244            psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
     245            pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL);
     246            psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
    244247        }
    245248
  • trunk/psphot/src/psphotEnsemblePSF.c

    r13006 r13035  
    9191
    9292        // fill in the model pixel values
     93        // XXX this function should just work with the (normalized) modelFlux pixels
     94        // XXX we should not have to create a complete copy
     95        // XXX review the use of the object mask: can we use a copied mask?
    9396        psImageInit (fitSource->pixels, 0.0);
    9497        psImageKeepCircle (fitSource->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
    95         pmModelAdd (fitSource->pixels, fitSource->mask, model, false, false);
     98        pmModelAdd (fitSource->pixels, fitSource->mask, model, PM_MODEL_OP_FULL);
    9699
    97100        // save source in list
     
    218221
    219222        // subtract object
    220         pmModelSub (Ri->pixels, Ri->mask, model, false, false);
     223        pmModelSub (Ri->pixels, Ri->mask, model, PM_MODEL_OP_FULL);
    221224        Ri->mode |= PM_SOURCE_MODE_SUBTRACTED;
    222225        if (!final) Ri->mode |= PM_SOURCE_MODE_TEMPSUB;
  • trunk/psphot/src/psphotFitSet.c

    r12792 r13035  
    11# include "psphotInternal.h"
    22
    3 // XXX this is not used in main psphot code
     3// This is not used in main psphot code (only in psphotModelTest.c)
    44bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode) {
    55
     
    2424    }
    2525
     26    // XXX pmSourceFitSet must cache the modelFlux?
    2627    pmSourceFitSet (source, modelSet, mode);
    2728
     
    3233    for (int i = 0; i < modelSet->n; i++) {
    3334        pmModel *model = modelSet->data[i];
    34         pmModelSub (source->pixels, source->mask, model, false, false);
     35        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
    3536
    3637        fprintf (stderr, "output parameters (obj %d):\n", i);
     
    4243    // write out
    4344    psphotSaveImage (NULL, source->pixels, "resid.fits");
    44     psphotSaveImage (NULL, source->mask, "mask.fits");
     45    psphotSaveImage (NULL, source->maskObj, "mask.fits");
    4546    return true;
    4647}
  • trunk/psphot/src/psphotGrowthCurve.c

    r12792 r13035  
    5353        radius = psf->growth->radius->data.F32[i];
    5454
     55        // NOTE: we use pmModelAdd not pmSourceAdd because we are not working with a normal source
    5556        psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK);
    56         pmModelAdd (image, mask, model, false, false);
     57        pmModelAdd (image, mask, model, PM_MODEL_OP_FULL);
    5758        pmSourcePhotometryAper (&apMag, model, image, mask);
    5859        psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK));
  • trunk/psphot/src/psphotGuessModels.c

    r12950 r13035  
    5656    source->modelPSF = modelPSF;
    5757    source->modelPSF->residuals = psf->residuals;
     58
     59    pmSourceCacheModel (source);
    5860  }
    5961  psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
  • trunk/psphot/src/psphotMakeResiduals.c

    r12950 r13035  
    7676        if (model == NULL) continue;  // model must be defined
    7777
    78         psImage *image  = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    79         psImage *weight = psImageCopy (NULL, source->weight, PS_TYPE_F32);
    80         psImage *mask   = psImageCopy (NULL, source->mask PS_TYPE_U8);
    81         pmModelSub (image, mask, model, false, false);
     78        psImage *image  = psImageCopy (NULL, source->pixels,   PS_TYPE_F32);
     79        psImage *weight = psImageCopy (NULL, source->weight,   PS_TYPE_F32);
     80        psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_U8);
     81        pmModelSub (image, mask, model, PM_MODEL_OP_FUNC);
    8282       
    8383        // re-normalize image and weight
     
    253253# endif
    254254
     255    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
     256
    255257    psFree (xC);
    256258    psFree (yC);
     
    272274    psphotSaveImage (NULL, resid->mask,   "resid.mk.fits");
    273275
    274     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
    275 
    276276    psf->residuals = resid;
    277277    return true;
  • trunk/psphot/src/psphotModelTest.c

    r12792 r13035  
    22static char DEFAULT_MODE[] = "EXT";
    33
     4// XXX consider this function : add more test information?
    45bool psphotModelTest (pmReadout *readout, psMetadata *recipe) {
    56
     
    158159
    159160    // define the pixels used for the fit
    160     psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PM_MASK_MARK);
     161    psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", PM_MASK_MARK);
    161162
    162163    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
     
    170171    // measure the source mags
    171172    pmSourcePhotometryModel (&fitMag, model);
    172     pmSourcePhotometryAper  (&obsMag, model, source->pixels, source->mask);
     173    pmSourcePhotometryAper  (&obsMag, model, source->pixels, source->maskObj);
    173174    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
    174175
     
    177178
    178179    // subtract object, leave local sky
    179     pmModelSub (source->pixels, source->mask, model, false, false);
     180    pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
    180181
    181182    fprintf (stderr, "output parameters: \n");
     
    186187    // write out
    187188    psphotSaveImage (NULL, source->pixels, "resid.fits");
    188     psphotSaveImage (NULL, source->mask, "mask.fits");
     189    psphotSaveImage (NULL, source->maskObj, "mask.fits");
    189190
    190191    exit (0);
  • trunk/psphot/src/psphotOutput.c

    r13008 r13035  
    5959        for (int j = 0; j < source->pixels->numCols; j++) {
    6060            // skip masked points
    61             if (source->mask->data.U8[i][j]) {
     61            if (source->maskObj->data.U8[i][j]) {
    6262                continue;
    6363            }
     
    7272                     source->pixels->data.F32[i][j],
    7373                     1.0 / source->weight->data.F32[i][j],
    74                      source->mask->data.U8[i][j]);
     74                     source->maskObj->data.U8[i][j]);
    7575        }
    7676    }
     
    104104    for (int i = 0; (sources != NULL) && (i < sources->n); i++) {
    105105        pmSource *source = (pmSource *) sources->data[i];
    106         pmModel *model = pmSourceSelectModel (source);
     106        pmModel *model = pmSourceGetModel (NULL, source);
    107107        if (model == NULL)
    108108            continue;
     
    154154    return header;
    155155}
    156 
    157 # if (0)
    158 // output functions: we have several fixed modes (sx, obj, cmp)
    159 void psphotOutput (pmReadout *readout, psMetadata *arguments) {
    160 
    161     bool status;
    162     char *outputFile = NULL;
    163 
    164     psMetadata *header = pmReadoutGetHeader (readout);
    165 
    166     pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    167     // sample PSF images??
    168     if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample);
    169     if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample);
    170 
    171     if (psfFile != NULL) {
    172         psMetadata *psfData = pmPSFtoMD (NULL, psf);
    173         psMetadataConfigWrite (psfData, psfFile);
    174         psFree (psfData);
    175     }
    176 }
    177 
    178 psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) {
    179 
    180     psRegion region = {x - dx, x + dx, y - dy, y + dy};
    181     psImage *view = psImageSubset (image, region);
    182     psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32);
    183     psImageInit (sample, 0);
    184     modelEXT->params->data.F32[2] = x;
    185     modelEXT->params->data.F32[3] = y;
    186     pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
    187     pmModelAdd (sample, NULL, modelPSF, false, false);
    188     psFree (modelPSF);
    189     return (sample);
    190 }
    191 
    192 bool psphotSamplePSFs (pmPSF *psf, psImage *image, char *output) {
    193 
    194     // make sample PSFs for 4 corners and the center
    195     psImage *sample;
    196 
    197     // optional dump of all rough source data
    198     if (output[0] == 0) return false;
    199 
    200     pmModel *modelEXT = pmModelAlloc (psf->type);
    201     modelEXT->params->data.F32[0] = 0;
    202     modelEXT->params->data.F32[1] = 1;
    203 
    204     psFits *fits = psFitsOpen (output, "w");
    205 
    206     // the centers are in parent coordinates; they do not need to correspond to valid pixels...
    207     sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25);
    208     psFitsWriteImage (fits, NULL, sample, 0);
    209     sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25);
    210     psFitsWriteImage (fits, NULL, sample, 0);
    211     sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25);
    212     psFitsWriteImage (fits, NULL, sample, 0);
    213     sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25);
    214     psFitsWriteImage (fits, NULL, sample, 0);
    215     sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25);
    216     psFitsWriteImage (fits, NULL, sample, 0);
    217 
    218     psFitsClose (fits);
    219     return (TRUE);
    220 }
    221 # endif
  • trunk/psphot/src/psphotRadialPlot.c

    r12900 r13035  
    6060    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    6161        for (int ix = 0; ix < source->pixels->numCols; ix++) {
    62             if (source->mask->data.U8[iy][ix]) {
     62            if (source->maskObj->data.U8[iy][ix]) {
    6363                rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
    6464                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
  • trunk/psphot/src/psphotRadiusChecks.c

    r12792 r13035  
    2525bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model)
    2626{
     27    psF32 *PAR = model->params->data.F32;
     28
     29    // XXX do we have a better value for the sky noise level?  not really...
    2730    pmMoments *moments = source->moments;
    28     // do we have a better value for the sky noise level?
    29     // not really...
    3031
    3132    // set the fit radius based on the object flux limit and the model
     
    4647    }
    4748
    48     bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit);
     49    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit);
     50
     51    // set the mask to flag the excluded pixels
     52    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", PM_MASK_MARK);
    4953    return status;
    5054}
    5155
    5256bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR) {
     57
     58    psF32 *PAR = model->params->data.F32;
    5359
    5460    pmMoments *moments = source->moments;
     
    6369    }
    6470
    65     bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit);
     71    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit);
     72
     73    // set the mask to flag the excluded pixels
     74    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", PM_MASK_MARK);
    6675    return status;
    6776}
     
    8695bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model) {
    8796
     97    psF32 *PAR = model->params->data.F32;
     98
    8899    pmMoments *moments = source->moments;
    89100    if (moments == NULL) return false;
     
    94105
    95106    // redefine the pixels if needed
    96     bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit);
     107    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit);
     108
     109    // set the mask to flag the excluded pixels
     110    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", PM_MASK_MARK);
    97111    return status;
    98112}
  • trunk/psphot/src/psphotReadout.c

    r13012 r13035  
    126126
    127127    // linear PSF fit to peaks
    128     psphotEnsemblePSF (readout, sources, recipe, psf, FALSE);
     128    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    129129    if (dump) psphotSaveImage (NULL, readout->image,  "image.v1.fits");
    130130    if (!strcasecmp (breakPt, "ENSEMBLE")) {
     
    141141
    142142    // linear PSF fit to remaining peaks
    143     psphotEnsemblePSF (readout, sources, recipe, psf, TRUE);
     143    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
    144144    if (dump) psphotSaveImage (NULL, readout->image,  "image.v4.fits");
    145145    if (!strcasecmp (breakPt, "PASS1")) {
     
    192192
    193193    // linear PSF fit to remaining peaks
    194     psphotEnsemblePSF (readout, sources, recipe, psf, TRUE);
     194    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
    195195    if (dump) psphotSaveImage (NULL, readout->image,  "image.v6.fits");
    196196
  • trunk/psphot/src/psphotReplaceUnfit.c

    r12792 r13035  
    11# include "psphotInternal.h"
    22
     3// replace the flux for sources which failed
    34bool psphotReplaceUnfit (psArray *sources) {
    45
     
    1516
    1617    replace:
    17         if (source->modelPSF == NULL) continue;
    18 
    19         psTrace ("psphot", 3, "replacing object at %f,%f\n",
    20                  source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS]);
    21 
    22         pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false);
     18        pmSourceAdd (source, PM_MODEL_OP_FULL);
    2319        source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
    2420        source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
     
    4036      if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
    4137
    42       // select appropriate model
    43       pmModel *model = pmSourceGetModel (NULL, source);
    44       if (model == NULL) continue;  // model must be defined
    45        
    46       psTrace ("psphot", 3, "replacing object at %f,%f\n",
    47                model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
    48 
    49       pmModelAdd (source->pixels, source->mask, model, false, false);
     38      pmSourceAdd (source, PM_MODEL_OP_FULL);
    5039      source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
    5140    }
     
    5443}
    5544
    56 // add or sub source replace or if the source has
     45// add source, if the source has been subtracted (or if we ignore the state)
    5746bool psphotAddWithTest (pmSource *source, bool useState) {
    5847
     
    6049    bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
    6150    if (state && useState) return true;
    62 
    63     // select appropriate model
    64     pmModel *model = pmSourceGetModel (NULL, source);
    65     if (model == NULL) return false;  // model must be defined
    66    
    67     psTrace ("psphot", 3, "replacing object at %f,%f\n",
    68              model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
    6951   
    7052    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
    7153    if (!state || !useState) {
    72         pmModelAdd (source->pixels, source->mask, model, false, false);
     54        pmSourceAdd (source, PM_MODEL_OP_FULL);
    7355    }
    7456    return true;
    7557}
    7658
    77 // add or sub source replace or if the source has
     59// sub source, if the source has been added (or if we ignore the state)
    7860bool psphotSubWithTest (pmSource *source, bool useState) {
    7961
     
    8264    if (state && useState) return true;
    8365
    84     // select appropriate model
    85     pmModel *model = pmSourceGetModel (NULL, source);
    86     if (model == NULL) return false;  // model must be defined
    87    
    88     psTrace ("psphot", 3, "replacing object at %f,%f\n",
    89              model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
    90    
    9166    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
    9267    if (!state || !useState) {
    93         pmModelSub (source->pixels, source->mask, model, false, false);
     68        pmSourceSub (source, PM_MODEL_OP_FULL);
    9469    }
    9570    return true;
     
    10277    bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
    10378    if (curState == newState) return true;
    104 
    105     // select appropriate model
    106     pmModel *model = pmSourceGetModel (NULL, source);
    107     if (model == NULL) return false;  // model must be defined
    108    
    109     psTrace ("psphot", 3, "replacing object at %f,%f\n",
    110              model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
    11179   
    11280    if (curState && !newState) {
    113         pmModelSub (source->pixels, source->mask, model, false, false);
     81        pmSourceSub (source, PM_MODEL_OP_FULL);
    11482    }
    11583    if (newState && !curState) {
    116         pmModelAdd (source->pixels, source->mask, model, false, false);
     84        pmSourceAdd (source, PM_MODEL_OP_FULL);
    11785    }
    11886    return true;
  • trunk/psphot/src/psphotSourceFits.c

    r12792 r13035  
    6767        // these should never be invalid values
    6868        // XXX drop these tests eventually
    69         if (isnan(model->params->data.F32[PM_PAR_I0])) psAbort("nan in blend fit");
     69        if (isnan(model->params->data.F32[PM_PAR_I0]))   psAbort("nan in blend fit");
    7070        if (isnan(model->params->data.F32[PM_PAR_XPOS])) psAbort("nan in blend fit");
    7171        if (isnan(model->params->data.F32[PM_PAR_YPOS])) psAbort("nan in blend fit");
     
    8383
    8484    // fit PSF model (set/unset the pixel mask)
    85     psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK);
    8685    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
    87     psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    8886
    8987    // correct model chisq for flux trend
     
    110108        blend->modelPSF = psMemIncrRefCounter (model);
    111109        psTrace ("psphot", 5, "fitted blend as PSF\n");
    112         pmModelSub (source->pixels, source->mask, model, false, false);
     110
     111        // build cached model and subtract
     112        pmSourceCacheModel (blend);
     113        pmSourceSub (blend, PM_MODEL_OP_FULL);
    113114        blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    114115        blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
     
    127128    psFree (sourceSet);
    128129
    129     psTrace ("psphot", 5, "fitted primary as PSF\n");
    130     pmModelSub (source->pixels, source->mask, PSF, false, false);
     130    // save the new, successful model
    131131    psFree (source->modelPSF);
    132132    source->modelPSF = PSF;
     133    psTrace ("psphot", 5, "fitted primary as PSF\n");
     134
     135    // build cached model and subtract
     136    pmSourceCacheModel (source);
     137    pmSourceSub (source, PM_MODEL_OP_FULL);
    133138    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    134139    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
     
    138143bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) {
    139144
    140     float x, y;
    141145    double chiTrend;
    142146
     
    150154    psphotCheckRadiusPSF (readout, source, PSF);
    151155
    152     x = PSF->params->data.F32[PM_PAR_XPOS];
    153     y = PSF->params->data.F32[PM_PAR_YPOS];
    154 
    155156    // fit PSF model (set/unset the pixel mask)
    156     psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK);
    157157    pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF);
    158     psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    159158
    160159    // correct model chisq for flux trend
     
    168167    }
    169168
    170     psTrace ("psphot", 5, "fitted as PSF\n");
    171     pmModelSub (source->pixels, source->mask, PSF, false, false);
    172 
    173169    // free old model, save new model
    174170    psFree (source->modelPSF);
    175171    source->modelPSF = PSF;
    176 
     172    psTrace ("psphot", 5, "fitted as PSF\n");
     173
     174    // build cached model and subtract
     175    pmSourceCacheModel (source);
     176    pmSourceSub (source, PM_MODEL_OP_FULL);
    177177    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    178178    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
     
    264264    // sub EXT
    265265    psFree (DBL);
    266     pmModelSub (source->pixels, source->mask, EXT, false, false);
    267     psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
    268266
    269267    // save new model
    270268    source->modelEXT = EXT;
     269
     270    // build cached model and subtract
     271    pmSourceCacheModel (source);
     272    pmSourceSub (source, PM_MODEL_OP_FULL);
     273    psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
     274
    271275    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    272276    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
     
    276280    // sub DLB
    277281    psFree (EXT);
    278     pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[0], false, false);
    279     pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[1], false, false);
    280     psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
    281282
    282283    // drop old model, save new second model...
     
    288289
    289290    // copy most data from the primary source (modelEXT, blends stay NULL)
    290     pmSource *newSrc = pmSourceAlloc ();
    291     newSrc->peak     = psMemIncrRefCounter (source->peak);
    292     newSrc->pixels   = psMemIncrRefCounter (source->pixels);
    293     newSrc->weight   = psMemIncrRefCounter (source->weight);
    294     newSrc->mask     = psMemIncrRefCounter (source->mask);
    295     newSrc->moments  = psMemIncrRefCounter (source->moments);
     291    // XXX use pmSourceCopy?
     292    pmSource *newSrc = pmSourceCopy (source);
    296293    newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
    297     newSrc->type     = source->type;
    298     newSrc->mode     = source->mode;
     294
     295    // build cached models and subtract
     296    pmSourceCacheModel (source);
     297    pmSourceSub (source, PM_MODEL_OP_FULL);
     298    pmSourceCacheModel (newSrc);
     299    pmSourceSub (newSrc, PM_MODEL_OP_FULL);
     300    psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
     301
    299302    psArrayAdd (sources, 100, newSrc);
    300303    psFree (newSrc);
     
    306309psArray *psphotFitDBL (pmReadout *readout, pmSource *source) {
    307310
    308     float x, y, dx, dy;
     311    float dx, dy;
    309312    pmModel *DBL;
    310313    pmModel *PSF;
     
    344347    modelSet->data[1] = DBL;
    345348
    346     x = source->moments->x;
    347     y = source->moments->y;
    348 
    349349    // fit PSF model (set/unset the pixel mask)
    350     psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK);
    351350    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
    352     psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    353 
    354351    return (modelSet);
    355352}
    356353
    357354pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) {
    358 
    359     float x, y;
    360355
    361356    NfitEXT ++;
     
    369364    psphotCheckRadiusEXT (readout, source, EXT);
    370365
    371     x = EXT->params->data.F32[PM_PAR_XPOS];
    372     y = EXT->params->data.F32[PM_PAR_YPOS];
    373 
    374366    if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) {
    375367        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy);
     
    377369
    378370    // fit EXT (not PSF) model (set/unset the pixel mask)
    379     psImageKeepCircle (source->mask, x, y, EXT->radiusFit, "OR", PM_MASK_MARK);
    380371    pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT);
    381     psImageKeepCircle (source->mask, x, y, EXT->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    382 
    383372    return (EXT);
    384373}
  • trunk/psphot/src/psphotTestPSF.c

    r13006 r13035  
    124124
    125125        // subtract model flux
    126         pmModelSub (source->pixels, source->mask, model, false, false);
     126        pmModelSub (source->pixels, source->mask, model, PM_MODEL_OP_FULL);
    127127    }
    128128    fclose (f);
  • trunk/psphot/src/psphotWeightBias.c

    r13006 r13035  
    55// only allow the normalization to vary
    66// XXX eventually, we should be able to do this with linear fitting...
     7// XXX UNUSED
    78bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
    89
     
    5960
    6061        // replace object in image
    61         pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false);
     62        pmModelAdd (source->pixels, source->mask, source->modelPSF, PM_MODEL_OP_FULL);
    6263
    6364        // make a temporary model (we don't keep the result of this analysis)
     
    7677
    7778        // re-subtract PSF for object, leave local sky
    78         pmModelSub (source->pixels, source->mask, source->modelPSF, false, false);
     79        pmModelSub (source->pixels, source->mask, source->modelPSF, PM_MODEL_OP_FULL);
    7980
    8081        PARp = source->modelPSF->params->data.F32;
Note: See TracChangeset for help on using the changeset viewer.