IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13900


Ignore:
Timestamp:
Jun 19, 2007, 4:40:44 PM (19 years ago)
Author:
Paul Price
Message:

Extensive changes to APIs to allow use of a nominated value to mask
against (the maskVal). Previously, the mask values were either
hard-coded (e.g., PM_MASK_SAT) or taken as anything non-zero. The
code is tested with psModules (which has similar changes) and does not
crash, but neither is it successful in marking all bad pixels. For
this reason, I have left the "gutter" pixels (cell gaps) set to 0
instead of NAN in pmFPAMosaic.

Location:
trunk/psphot/src
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphot.h

    r13593 r13900  
    1616psString        psphotVersionLong(void);
    1717
    18 bool            psphotModelTest (pmReadout *readout, psMetadata *recipe);
     18bool            psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark);
    1919bool            psphotReadout (pmConfig *config, const pmFPAview *view);
    2020bool            psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources);
    … …  
    2626
    2727// psphotReadout functions
    28 bool            psphotImageMedian (pmConfig *config, const pmFPAview *view);
     28bool            psphotImageMedian (pmConfig *config, const pmFPAview *view, psMaskType maskVal);
    2929psArray        *psphotFindPeaks (pmReadout *readout, psMetadata *recipe,
    30                                  const bool returnFootprints, const int pass);
     30                                 const bool returnFootprints, const int pass, psMaskType maskVal);
    3131#include "pmFootprint.h"
    3232psErrorCode     psphotCullPeaks(const psImage *img, const psImage *weight,
    3333                                const psMetadata *recipe, psArray *footprints);
    34 psArray        *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *allpeaks);
    35 bool            psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump);
     34psArray        *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *allpeaks, psMaskType maskVal, psMaskType mark);
     35bool            psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump, psMaskType maskSat);
    3636bool            psphotBasicDeblend (psArray *sources, psMetadata *recipe);
    37 pmPSF          *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe);
     37pmPSF          *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark);
    3838bool            psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
    3939bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
    … …  
    4141bool            psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
    4242#endif
    43 bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
    44 bool            psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
    45 bool            psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
    46 bool            psphotReplaceUnfit (psArray *sources);
    47 bool            psphotReplaceAll (psArray *sources);
    48 bool            psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
    49 bool            psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background);
     43bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final, psMaskType maskVal);
     44bool            psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal);
     45bool            psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal);
     46bool            psphotReplaceUnfit (psArray *sources, psMaskType maskVal);
     47bool            psphotReplaceAll (psArray *sources, psMaskType maskVal);
     48bool            psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark);
     49bool            psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background, psMaskType maskVal, psMaskType mark);
    5050bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view);
    5151
    … …  
    5656int             pmSourceSortBySN (const void **a, const void **b);
    5757int             pmSourceSortByY (const void **a, const void **b);
    58 bool            psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore);
    59 bool            psphotMaskReadout (pmReadout *readout, psMetadata *recipe, pmConfig *config);
     58bool            psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psMaskType maskVal);
     59bool            psphotMaskReadout (pmReadout *readout, psMetadata *recipe, psMaskType maskVal);
    6060void            psphotSourceFreePixels (psArray *sources);
    6161
    … …  
    8383bool            psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout);
    8484bool            psphotInitLimitsEXT (psMetadata *recipe);
    85 bool            psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf);
    86 bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf);
    87 bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf);
    88 pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source);
    89 psArray        *psphotFitDBL (pmReadout *readout, pmSource *source);
     85bool            psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal);
     86bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal);
     87bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal);
     88pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal);
     89psArray        *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal);
    9090
    9191// functions to support simultaneous multi-source fitting
    92 bool            psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode);
     92bool            psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode, psMaskType maskVal);
    9393
    9494// plotting functions (available if libkapa is installed)
    … …  
    101101pmPSF          *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe);
    102102bool            psphotSetHeaderNstars (psMetadata *recipe, psArray *sources);
    103 bool            psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add);
     103bool            psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add, psMaskType maskVal);
    104104bool            psphotRadialPlot (int *kapa, const char *filename, pmSource *source);
    105 bool            psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe);
     105bool            psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal);
    106106bool            psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY);
    107107
    108 bool            psphotAddWithTest (pmSource *source, bool useState);
    109 bool            psphotSubWithTest (pmSource *source, bool useState);
    110 bool            psphotSetState (pmSource *source, bool curState);
     108bool            psphotAddWithTest (pmSource *source, bool useState, psMaskType maskVal);
     109bool            psphotSubWithTest (pmSource *source, bool useState, psMaskType maskVal);
     110bool            psphotSetState (pmSource *source, bool curState, psMaskType maskVal);
    111111bool            psphotDeblendSatstars (psArray *sources, psMetadata *recipe);
    112112bool            psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe);
    113113
    114 bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf);
     114bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal);
    115115
    116116#endif

  • trunk/psphot/src/psphotAddNoise.c

    r13330 r13900  
    11# include "psphotInternal.h"
    22
    3 bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) {
     3bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add, psMaskType maskVal) {
    44
    55    bool status = false;
    … …  
    3434    // loop over all source
    3535    for (int i = 0; i < sources->n; i++) {
    36         pmSource *source = sources->data[i];
     36        pmSource *source = sources->data[i];
    3737
    38         // skip sources which were not subtracted
    39         if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
     38        // skip sources which were not subtracted
     39        if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
    4040
    41         // select appropriate model
    42         pmModel *model = pmSourceGetModel (NULL, source);
    43         if (model == NULL) continue;  // model must be defined
    44        
    45         if (add) {
    46             psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
    47         } else {
    48             psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
    49         }
     41        // select appropriate model
     42        pmModel *model = pmSourceGetModel (NULL, source);
     43        if (model == NULL) continue;  // model must be defined
    5044
    51         psF32 *PAR = model->params->data.F32;
     45        if (add) {
     46            psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     47        } else {
     48            psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     49        }
    5250
    53         // save original values
    54         float oldI0  = PAR[PM_PAR_I0];
    55         oldshape.sx  = PAR[PM_PAR_SXX];
    56         oldshape.sy  = PAR[PM_PAR_SYY];
    57         oldshape.sxy = PAR[PM_PAR_SXY];
     51        psF32 *PAR = model->params->data.F32;
    5852
    59         // increase size and height of source
    60         axes = psEllipseShapeToAxes (oldshape, 20.0);
    61         axes.major *= SIZE;
    62         axes.minor *= SIZE;
    63         newshape = psEllipseAxesToShape (axes);
    64         PAR[PM_PAR_I0]  = FACTOR*PS_SQR(oldI0);
    65         PAR[PM_PAR_SXX] = newshape.sx;
    66         PAR[PM_PAR_SYY] = newshape.sy;
    67         PAR[PM_PAR_SXY] = newshape.sxy;
     53        // save original values
     54        float oldI0  = PAR[PM_PAR_I0];
     55        oldshape.sx  = PAR[PM_PAR_SXX];
     56        oldshape.sy  = PAR[PM_PAR_SYY];
     57        oldshape.sxy = PAR[PM_PAR_SXY];
    6858
    69         // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected
    70         pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add);
    71        
    72         // restore original values
    73         PAR[PM_PAR_I0]  = oldI0;       
    74         PAR[PM_PAR_SXX] = oldshape.sx;
    75         PAR[PM_PAR_SYY] = oldshape.sy;
    76         PAR[PM_PAR_SXY] = oldshape.sxy;
     59        // increase size and height of source
     60        axes = psEllipseShapeToAxes (oldshape, 20.0);
     61        axes.major *= SIZE;
     62        axes.minor *= SIZE;
     63        newshape = psEllipseAxesToShape (axes);
     64        PAR[PM_PAR_I0]  = FACTOR*PS_SQR(oldI0);
     65        PAR[PM_PAR_SXX] = newshape.sx;
     66        PAR[PM_PAR_SYY] = newshape.sy;
     67        PAR[PM_PAR_SXY] = newshape.sxy;
     68
     69        // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected
     70        pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add, maskVal);
     71
     72        // restore original values
     73        PAR[PM_PAR_I0]  = oldI0;
     74        PAR[PM_PAR_SXX] = oldshape.sx;
     75        PAR[PM_PAR_SYY] = oldshape.sy;
     76        PAR[PM_PAR_SXY] = oldshape.sxy;
    7777    }
    7878    if (add) {
    79         psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
     79        psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
    8080    } else {
    81         psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
     81        psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
    8282    }
    8383    return true;

  • trunk/psphot/src/psphotApResid.c

    r13864 r13900  
    66
    77// measure the aperture residual statistics
    8 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf)
     8bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark)
    99{
    1010    int Nfail = 0;
    … …  
    3131    bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
    3232    int NSTAR_APERTURE_CORRECTION_MIN =
    33         psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN");
     33        psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN");
    3434    if (!status) {
    35         NSTAR_APERTURE_CORRECTION_MIN = 5;
     35        NSTAR_APERTURE_CORRECTION_MIN = 5;
    3636    }
    3737
    … …  
    4747    psf->growth = pmGrowthCurveAlloc (PSF_FIT_PAD, 100.0, REF_RADIUS);
    4848
    49     if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH)) {
    50         psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
    51         psFree(psf->growth); psf->growth = NULL;
    52         return false;
     49    if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, mark)) {
     50        psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
     51        psFree(psf->growth); psf->growth = NULL;
     52        return false;
    5353    }
    5454
    … …  
    8080        // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
    8181        // will fail if below S/N threshold or model is missing
    82         if (!pmSourceMagnitudes (source, psf, photMode)) {
     82        if (!pmSourceMagnitudes (source, psf, photMode, maskVal, mark)) {
    8383            Nskip ++;
    8484            continue;
    … …  
    135135    if (Npsf < NSTAR_APERTURE_CORRECTION_MIN) {
    136136        psError(PSPHOT_ERR_APERTURE, true, "Only %d valid aperture residual sources (need %d), giving up",
    137                 Npsf, NSTAR_APERTURE_CORRECTION_MIN);
     137                Npsf, NSTAR_APERTURE_CORRECTION_MIN);
    138138        return false;
    139139    }
    … …  
    154154    stats->max = 3.0;
    155155
    156 #define P_APTREND_SWITCH_CLEANUP        /* Cleanup memory on error in ApTrendOption switch */ \
     156#define P_APTREND_SWITCH_CLEANUP        /* Cleanup memory on error in ApTrendOption switch */ \
    157157    psFree(psf->growth); psf->growth = NULL; \
    158158    psFree(mask); \
    … …  
    172172        break;
    173173      case PM_PSF_APTREND_CONSTANT:
    174         stats->clipIter = 2;
    175         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    176         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    177             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    178             P_APTREND_SWITCH_CLEANUP;
    179             return false;
    180         }
    181         // apply the fit
    182         stats->clipIter = 3;
    183         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    184         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    185             psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
    186             P_APTREND_SWITCH_CLEANUP;
    187             return false;
    188         }
    189         break;
     174        stats->clipIter = 2;
     175        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     176        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     177            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     178            P_APTREND_SWITCH_CLEANUP;
     179            return false;
     180        }
     181        // apply the fit
     182        stats->clipIter = 3;
     183        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     184        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     185            psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
     186            P_APTREND_SWITCH_CLEANUP;
     187            return false;
     188        }
     189        break;
    190190      case PM_PSF_APTREND_SKYBIAS:
    191         // first clip out objects which are too far from the median
    192         stats->clipIter = 2;
    193         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    194         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    195             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    196             P_APTREND_SWITCH_CLEANUP;
    197             return false;
    198         }
    199         // apply the fit
    200         stats->clipIter = 3;
    201         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    202         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    203             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    204             P_APTREND_SWITCH_CLEANUP;
    205             return false;
    206         }
    207         break;
     191        // first clip out objects which are too far from the median
     192        stats->clipIter = 2;
     193        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     194        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     195            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     196            P_APTREND_SWITCH_CLEANUP;
     197            return false;
     198        }
     199        // apply the fit
     200        stats->clipIter = 3;
     201        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     202        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     203            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     204            P_APTREND_SWITCH_CLEANUP;
     205            return false;
     206        }
     207        break;
    208208      case PM_PSF_APTREND_SKYSAT:
    209         // first clip out objects which are too far from the median
    210         stats->clipIter = 2;
    211         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    212         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    213             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    214             P_APTREND_SWITCH_CLEANUP;
    215             return false;
    216         }
    217         // apply the fit
    218         stats->clipIter = 2;
    219         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    220         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    221             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    222             P_APTREND_SWITCH_CLEANUP;
    223             return false;
    224         }
    225         // apply the fit
    226         stats->clipIter = 3;
    227         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
    228         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    229             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
    230             P_APTREND_SWITCH_CLEANUP;
    231             return false;
    232         }
    233         break;
     209        // first clip out objects which are too far from the median
     210        stats->clipIter = 2;
     211        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     212        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     213            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     214            P_APTREND_SWITCH_CLEANUP;
     215            return false;
     216        }
     217        // apply the fit
     218        stats->clipIter = 2;
     219        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     220        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     221            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     222            P_APTREND_SWITCH_CLEANUP;
     223            return false;
     224        }
     225        // apply the fit
     226        stats->clipIter = 3;
     227        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
     228        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     229            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
     230            P_APTREND_SWITCH_CLEANUP;
     231            return false;
     232        }
     233        break;
    234234      case PM_PSF_APTREND_XY_LIN:
    235         // first clip out objects which are too far from the median
    236         stats->clipIter = 2;
    237         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    238         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    239             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    240             P_APTREND_SWITCH_CLEANUP;
    241             return false;
    242         }
    243         // apply the fit
    244         stats->clipIter = 3;
    245         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
    246         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    247             psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
    248             P_APTREND_SWITCH_CLEANUP;
    249             return false;
    250         }
    251         break;
     235        // first clip out objects which are too far from the median
     236        stats->clipIter = 2;
     237        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     238        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     239            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     240            P_APTREND_SWITCH_CLEANUP;
     241            return false;
     242        }
     243        // apply the fit
     244        stats->clipIter = 3;
     245        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
     246        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     247            psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
     248            P_APTREND_SWITCH_CLEANUP;
     249            return false;
     250        }
     251        break;
    252252      case PM_PSF_APTREND_XY_QUAD:
    253         // first clip out objects which are too far from the median
    254         stats->clipIter = 2;
    255         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    256         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    257             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    258             P_APTREND_SWITCH_CLEANUP;
    259             return false;
    260         }
    261         // apply the fit
    262         stats->clipIter = 3;
    263         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
    264         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    265             psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
    266             P_APTREND_SWITCH_CLEANUP;
    267             return false;
    268         }
    269         break;
     253        // first clip out objects which are too far from the median
     254        stats->clipIter = 2;
     255        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     256        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     257            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     258            P_APTREND_SWITCH_CLEANUP;
     259            return false;
     260        }
     261        // apply the fit
     262        stats->clipIter = 3;
     263        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
     264        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     265            psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
     266            P_APTREND_SWITCH_CLEANUP;
     267            return false;
     268        }
     269        break;
    270270      case PM_PSF_APTREND_SKY_XY_LIN:
    271         // first clip out objects which are too far from the median
    272         stats->clipIter = 2;
    273         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    274         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    275             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    276             P_APTREND_SWITCH_CLEANUP;
    277             return false;
    278         }
    279         // apply the fit
    280         stats->clipIter = 3;
    281         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
    282         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    283             psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
    284             P_APTREND_SWITCH_CLEANUP;
    285             return false;
    286         }
    287         break;
     271        // first clip out objects which are too far from the median
     272        stats->clipIter = 2;
     273        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     274        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     275            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     276            P_APTREND_SWITCH_CLEANUP;
     277            return false;
     278        }
     279        // apply the fit
     280        stats->clipIter = 3;
     281        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
     282        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     283            psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
     284            P_APTREND_SWITCH_CLEANUP;
     285            return false;
     286        }
     287        break;
    288288      case PM_PSF_APTREND_SKYSAT_XY_LIN:
    289         // first clip out objects which are too far from the median
    290         stats->clipIter = 2;
    291         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    292         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    293             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    294             P_APTREND_SWITCH_CLEANUP;
    295             return false;
    296         }
    297         // apply the fit
    298         stats->clipIter = 3;
    299         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    300         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    301             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    302             P_APTREND_SWITCH_CLEANUP;
    303             return false;
    304         }
    305         // apply the fit
    306         stats->clipIter = 3;
    307         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
    308         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    309             psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
    310             P_APTREND_SWITCH_CLEANUP;
    311             return false;
    312         }
    313         break;
     289        // first clip out objects which are too far from the median
     290        stats->clipIter = 2;
     291        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     292        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     293            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     294            P_APTREND_SWITCH_CLEANUP;
     295            return false;
     296        }
     297        // apply the fit
     298        stats->clipIter = 3;
     299        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     300        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     301            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     302            P_APTREND_SWITCH_CLEANUP;
     303            return false;
     304        }
     305        // apply the fit
     306        stats->clipIter = 3;
     307        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
     308        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     309            psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
     310            P_APTREND_SWITCH_CLEANUP;
     311            return false;
     312        }
     313        break;
    314314      case PM_PSF_APTREND_ALL:
    315         // first clip out objects which are too far from the median
    316         stats->clipIter = 2;
    317         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    318         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    319             psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
    320             P_APTREND_SWITCH_CLEANUP;
    321             return false;
    322         }
    323         // fit just SkyBias and clip out objects which are too far from the median
    324         stats->clipIter = 2;
    325         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    326         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    327             psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
    328             P_APTREND_SWITCH_CLEANUP;
    329             return false;
    330         }
    331         // finally, fit x, y, SkyBias and clip out objects which are too far from the median
    332         stats->clipIter = 3;
    333         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
    334         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    335             psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
    336             P_APTREND_SWITCH_CLEANUP;
    337             return false;
    338         }
    339         break;
     315        // first clip out objects which are too far from the median
     316        stats->clipIter = 2;
     317        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     318        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     319            psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
     320            P_APTREND_SWITCH_CLEANUP;
     321            return false;
     322        }
     323        // fit just SkyBias and clip out objects which are too far from the median
     324        stats->clipIter = 2;
     325        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     326        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     327            psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
     328            P_APTREND_SWITCH_CLEANUP;
     329            return false;
     330        }
     331        // finally, fit x, y, SkyBias and clip out objects which are too far from the median
     332        stats->clipIter = 3;
     333        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
     334        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
     335            psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
     336            P_APTREND_SWITCH_CLEANUP;
     337            return false;
     338        }
     339        break;
    340340      default:
    341341        psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
    … …  
    343343    }
    344344#undef P_APTREND_SWITCH_CLEANUP
    345    
     345
    346346    // construct the fitted values and the residuals
    347347    psVector *fit   = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux);
    … …  
    376376    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
    377377
    378     psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n", 
    379               Nkeep, Npsf, psTimerMark ("psphot"));
     378    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n",
     379              Nkeep, Npsf, psTimerMark ("psphot"));
    380380    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f : %f bias, %f skysat\n",
    381381              psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat);

  • trunk/psphot/src/psphotBlendFit.c

    r13035 r13900  
    22
    33// XXX I don't like this name
    4 bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
     4bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) {
    55
    66    int Nfit = 0;
    … …  
    1414    // source analysis is done in S/N order (brightest first)
    1515    sources = psArraySort (sources, pmSourceSortBySN);
    16    
     16
    1717    // S/N limit to perform full non-linear fits
    1818    float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
    … …  
    3030
    3131    for (int i = 0; i < sources->n; i++) {
    32         // if (i%100 == 0) psphotFitSummary ();
     32        // if (i%100 == 0) psphotFitSummary ();
    3333
    34         pmSource *source = sources->data[i];
     34        pmSource *source = sources->data[i];
    3535
    36         // skip non-astronomical objects (very likely defects)
    37         if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
    38         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    39         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     36        // skip non-astronomical objects (very likely defects)
     37        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
     38        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     39        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    4040
    41         // skip DBL second sources (ie, added by psphotFitBlob)
    42         if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
     41        // skip DBL second sources (ie, added by psphotFitBlob)
     42        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
    4343
    44         // limit selection to some SN limit
    45         // XXX this should use peak?
    46         if (source->moments == NULL) continue;
    47         if (source->moments->SN < FIT_SN_LIM) continue;
     44        // limit selection to some SN limit
     45        // XXX this should use peak?
     46        if (source->moments == NULL) continue;
     47        if (source->moments->SN < FIT_SN_LIM) continue;
    4848
    49         // XXX this should use peak?
    50         if (source->moments->x < AnalysisRegion.x0) continue;
    51         if (source->moments->y < AnalysisRegion.y0) continue;
    52         if (source->moments->x > AnalysisRegion.x1) continue;
    53         if (source->moments->y > AnalysisRegion.y1) continue;
     49        // XXX this should use peak?
     50        if (source->moments->x < AnalysisRegion.x0) continue;
     51        if (source->moments->y < AnalysisRegion.y0) continue;
     52        if (source->moments->x > AnalysisRegion.x1) continue;
     53        if (source->moments->y > AnalysisRegion.y1) continue;
    5454
    55         // if model is NULL, we don't have a starting guess
    56         if (source->modelPSF == NULL) continue;
     55        // if model is NULL, we don't have a starting guess
     56        if (source->modelPSF == NULL) continue;
    5757
    58         // skip sources which are insignificant flux?
    59         if (source->modelPSF->params->data.F32[1] < 0.1) {
    60             psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
    61                      source->modelPSF->params->data.F32[1],
    62                      source->modelPSF->params->data.F32[2],
    63                      source->modelPSF->params->data.F32[3]);
    64             continue;
    65         }
     58        // skip sources which are insignificant flux?
     59        if (source->modelPSF->params->data.F32[1] < 0.1) {
     60            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
     61                     source->modelPSF->params->data.F32[1],
     62                     source->modelPSF->params->data.F32[2],
     63                     source->modelPSF->params->data.F32[3]);
     64            continue;
     65        }
    6666
    67         // replace object in image
    68         if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
    69             pmSourceAdd (source, PM_MODEL_OP_FULL);
    70         }
    71         Nfit ++;
     67        // replace object in image
     68        if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
     69            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     70        }
     71        Nfit ++;
    7272
    73         // try fitting PSFs, then try extended sources
    74         // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
    75         if (psphotFitBlend (readout, source, psf)) {
    76             psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y);
    77             Npsf ++;
    78             continue;
    79         }
    80         if (psphotFitBlob (readout, source, sources, psf)) {
    81             psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y);
    82             Next ++;
    83             continue;
    84         }
     73        // try fitting PSFs, then try extended sources
     74        // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
     75        if (psphotFitBlend (readout, source, psf, maskVal)) {
     76            psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y);
     77            Npsf ++;
     78            continue;
     79        }
     80        if (psphotFitBlob (readout, source, sources, psf, maskVal)) {
     81            psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y);
     82            Next ++;
     83            continue;
     84        }
    8585
    86         psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y);
    87         Nfail ++;
     86        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y);
     87        Nfail ++;
    8888
    89         // re-subtract the object, leave local sky
    90         pmSourceCacheModel (source);
    91         pmSourceSub (source, PM_MODEL_OP_FULL);
    92         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
    93         source->mode |= PM_SOURCE_MODE_TEMPSUB;
     89        // re-subtract the object, leave local sky
     90        pmSourceCacheModel (source, maskVal);
     91        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     92        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     93        source->mode |= PM_SOURCE_MODE_TEMPSUB;
    9494    }
    9595

  • trunk/psphot/src/psphotChoosePSF.c

    r13834 r13900  
    22
    33// try PSF models and select best option
    4 pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) {
     4pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {
    55
    66    bool status;
    … …  
    4343    psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md);
    4444    if (!psfTrendMask) {
    45         psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
    46         return NULL;
     45        psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
     46        return NULL;
    4747    }
    4848
    … …  
    5353        pmSource *source = sources->data[i];
    5454        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
    55             // keep NSTARS PSF stars, unmark the rest
    56             if (stars->n < NSTARS) {
    57                 psArrayAdd (stars, 200, source);
    58             } else {
    59                 source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
    60             }
    61         }
     55            // keep NSTARS PSF stars, unmark the rest
     56            if (stars->n < NSTARS) {
     57                psArrayAdd (stars, 200, source);
     58            } else {
     59                source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
     60            }
     61        }
    6262    }
    6363    psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);
    … …  
    6565    if (stars->n == 0) {
    6666        psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates");
    67         psFree (stars);
    68         psFree (psfTrendMask);
     67        psFree (stars);
     68        psFree (psfTrendMask);
    6969        return NULL;
    7070    }
    … …  
    9191        psMetadataItem *item = psListGetAndIncrement (iter);
    9292        char *modelName = item->data.V;
    93         models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS);
     93        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark);
    9494    }
    9595
    … …  
    127127    // print/dump psf parameters
    128128    if (psTraceGetLevel("psphot") >= 5) {
    129         for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) {
    130             psPolynomial2D *poly = try->psf->params_NEW->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             }
    137         }
     129        for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) {
     130            psPolynomial2D *poly = try->psf->params_NEW->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            }
     137        }
    138138    }
    139139
    … …  
    147147    psVector *dSN = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
    148148    for (int i = 0; i < try->sources->n; i++) {
    149         // masked for: bad model fit, outlier in parameters
    150         if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
    151             continue;
    152 
    153         pmSource *source = try->sources->data[i];
    154         Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX];
    155         Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY];
    156         dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0];
    157         dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX];
    158         dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY];
    159         Sx->n ++;
    160         Sy->n ++;
    161         dSN->n ++;
    162         dSx->n ++;
    163         dSy->n ++;
     149        // masked for: bad model fit, outlier in parameters
     150        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
     151            continue;
     152
     153        pmSource *source = try->sources->data[i];
     154        Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX];
     155        Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY];
     156        dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0];
     157        dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX];
     158        dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY];
     159        Sx->n ++;
     160        Sy->n ++;
     161        dSN->n ++;
     162        dSx->n ++;
     163        dSy->n ++;
    164164    }
    165165    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    … …  
    174174
    175175    for (int i = 0; i < dSN->n; i++) {
    176         dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]);
    177         dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]);
     176        dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]);
     177        dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]);
    178178    }
    179179
    … …  
    206206
    207207    // build a PSF residual image
    208     if (!psphotMakeResiduals (try->sources, recipe, try->psf)) {
    209         psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF");
    210         psFree (models);
    211         return NULL;
     208    if (!psphotMakeResiduals (try->sources, recipe, try->psf, maskVal)) {
     209        psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF");
     210        psFree (models);
     211        return NULL;
    212212    }
    213213
    214214    // XXX test dump of psf star data and psf-subtracted image
    215     if (psTraceGetLevel("psphot.psfstars") > 5) { 
    216         psphotSaveImage (NULL, readout->image,  "rawstars.fits");
    217 
    218         for (int i = 0; i < try->sources->n; i++) {
    219             // masked for: bad model fit, outlier in parameters
    220             if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
    221                 continue;
    222 
    223             pmSource *source = try->sources->data[i];
    224             float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    225             float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    226 
    227             // set the mask and subtract the PSF model
    228             // XXX should we be using maskObj? should we be unsetting the mask?
    229             // use pmModelSub because modelFlux has not been generated
    230             assert (source->maskObj);
    231             psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
    232             pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL);
    233             psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
    234         }
    235 
    236         FILE *f = fopen ("shapes.dat", "w");
    237         for (int i = 0; i < try->sources->n; i++) {
    238             psF32 inPar[10];
    239 
    240             // masked for: bad model fit, outlier in parameters
    241             if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    242 
    243             pmSource *source = try->sources->data[i];
    244             psF32 *outPar = source->modelEXT->params->data.F32;
    245 
    246             psEllipseShape shape;
    247 
    248             shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
    249             shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
    250             shape.sxy = outPar[PM_PAR_SXY];
    251 
    252             psEllipsePol pol = pmPSF_ModelToFit (outPar);
    253             inPar[PM_PAR_E0] = pol.e0;
    254             inPar[PM_PAR_E1] = pol.e1;
    255             inPar[PM_PAR_E2] = pol.e2;
    256             pmPSF_FitToModel (inPar, 0.1);
    257 
    258             psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
    259             psEllipseAxes axes2;
    260             (void)psEllipsePolToAxes(pol, 0.1, &axes2);
    261             psEllipsePol pol2 = psEllipseAxesToPol (axes1);
    262 
    263             fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
    264                      i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
    265                      outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
    266                      pol.e0, pol.e1, pol.e2,
    267                      pol2.e0, pol2.e1, pol2.e2,
    268                      inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
    269                      axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
    270                      axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
    271                      );
    272         }
    273         fclose (f);
    274 
    275         psphotSaveImage (NULL, readout->image,  "psfstars.fits");
    276         pmSourcesWritePSFs (try->sources, "psfstars.dat");
    277         pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
    278         psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
    279         psMetadataConfigWrite (psfData, "psfmodel.dat");
    280         psFree (psfData);
    281         psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
    282         exit (0);
     215    if (psTraceGetLevel("psphot.psfstars") > 5) {
     216        psphotSaveImage (NULL, readout->image,  "rawstars.fits");
     217
     218        for (int i = 0; i < try->sources->n; i++) {
     219            // masked for: bad model fit, outlier in parameters
     220            if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
     221                continue;
     222
     223            pmSource *source = try->sources->data[i];
     224            float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     225            float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     226
     227            // set the mask and subtract the PSF model
     228            // XXX should we be using maskObj? should we be unsetting the mask?
     229            // use pmModelSub because modelFlux has not been generated
     230            assert (source->maskObj);
     231            psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
     232            pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
     233            psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
     234        }
     235
     236        FILE *f = fopen ("shapes.dat", "w");
     237        for (int i = 0; i < try->sources->n; i++) {
     238            psF32 inPar[10];
     239
     240            // masked for: bad model fit, outlier in parameters
     241            if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
     242
     243            pmSource *source = try->sources->data[i];
     244            psF32 *outPar = source->modelEXT->params->data.F32;
     245
     246            psEllipseShape shape;
     247
     248            shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
     249            shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
     250            shape.sxy = outPar[PM_PAR_SXY];
     251
     252            psEllipsePol pol = pmPSF_ModelToFit (outPar);
     253            inPar[PM_PAR_E0] = pol.e0;
     254            inPar[PM_PAR_E1] = pol.e1;
     255            inPar[PM_PAR_E2] = pol.e2;
     256            pmPSF_FitToModel (inPar, 0.1);
     257
     258            psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
     259            psEllipseAxes axes2;
     260            (void)psEllipsePolToAxes(pol, 0.1, &axes2);
     261            psEllipsePol pol2 = psEllipseAxesToPol (axes1);
     262
     263            fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
     264                     i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
     265                     outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
     266                     pol.e0, pol.e1, pol.e2,
     267                     pol2.e0, pol2.e1, pol2.e2,
     268                     inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
     269                     axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
     270                     axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
     271                     );
     272        }
     273        fclose (f);
     274
     275        psphotSaveImage (NULL, readout->image,  "psfstars.fits");
     276        pmSourcesWritePSFs (try->sources, "psfstars.dat");
     277        pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
     278        psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
     279        psMetadataConfigWrite (psfData, "psfmodel.dat");
     280        psFree (psfData);
     281        psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
     282        exit (0);
    283283    }
    284284
    … …  
    290290
    291291    if (!psphotPSFstats (readout, recipe, psf)) {
    292         psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
    293         psFree(psf);
    294         return NULL;
     292        psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
     293        psFree(psf);
     294        return NULL;
    295295    }
    296296
    … …  
    327327    pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
    328328    if (modelPSF == NULL) {
    329         psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre");
    330         psFree(modelEXT);
    331         return false;
     329        psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre");
     330        psFree(modelEXT);
     331        return false;
    332332    }
    333333
    … …  
    346346    psF64 FWHM_Y = FWHM_X * (axes.minor / axes.major);
    347347
    348     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X",   PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
    349     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y",   PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
    350     psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
     348    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X",   PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
     349    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y",   PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
     350    psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
    351351    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
    352352    psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
    … …  
    386386        moments.xy = source->moments->Sxy;
    387387
    388         // limit axis ratio < 20.0
     388        // limit axis ratio < 20.0
    389389        axes = psEllipseMomentsToAxes (moments, 20.0);
    390390

  • trunk/psphot/src/psphotFindPeaks.c

    r13430 r13900  
    11# include "psphotInternal.h"
    22
    3 // In this function, we smooth the image, then search for the peaks 
     3// In this function, we smooth the image, then search for the peaks
    44psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe,
    5                           bool returnFootprints,
    6                           const int pass) {
     5                          bool returnFootprints,
     6                          const int pass, psMaskType maskVal) {
    77
    88    float SIGMA_SMTH, NSIGMA_SMTH, NSIGMA_PEAK;
    … …  
    1313
    1414    if (pass == 1) {
    15         SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
    16         NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
     15        SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
     16        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    1717    } else {
    18         bool status_x, status_y;
    19         float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
    20         float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
    21         if (!status_x | !status_y) {
    22             psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined");
    23             return false;
    24         }
    25         SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
    26         NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
     18        bool status_x, status_y;
     19        float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
     20        float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
     21        if (!status_x | !status_y) {
     22            psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined");
     23            return false;
     24        }
     25        SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
     26        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    2727    }
    2828
    2929    // smooth the image, applying the mask as we go
    3030    psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
    31     psImageSmoothMaskF32 (smooth_im, readout->mask, 0xff, SIGMA_SMTH, NSIGMA_SMTH);
     31    psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH);
    3232    psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot"));
    3333
    3434    // smooth the weight, applying the mask as we go
    3535    psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
    36     psImageSmoothMaskF32 (smooth_wt, readout->mask, 0xff, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
     36    psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
    3737    psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot"));
    3838
    3939    psImage *mask = readout->mask;
    4040
    41     // optionally save example images under trace 
     41    // optionally save example images under trace
    4242    if (psTraceGetLevel("psphot") > 5) {
    43         char name[64];
    44         sprintf (name, "imsmooth.v%d.fits", pass);
    45         psphotSaveImage (NULL, smooth_im, name);
    46         sprintf (name, "wtsmooth.v%d.fits", pass);
    47         psphotSaveImage (NULL, smooth_wt, name);
     43        char name[64];
     44        sprintf (name, "imsmooth.v%d.fits", pass);
     45        psphotSaveImage (NULL, smooth_im, name);
     46        sprintf (name, "wtsmooth.v%d.fits", pass);
     47        psphotSaveImage (NULL, smooth_wt, name);
    4848    }
    4949
    5050    // build the significance image on top of smooth_im
    5151    for (int j = 0; j < smooth_im->numRows; j++) {
    52         for (int i = 0; i < smooth_im->numCols; i++) {
    53             float value = smooth_im->data.F32[j][i];
    54             if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || mask->data.U8[j][i]) {
    55                 smooth_im->data.F32[j][i] = 0.0;
    56             } else {
    57                 smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i];
    58             }
    59         }
     52        for (int i = 0; i < smooth_im->numCols; i++) {
     53            float value = smooth_im->data.F32[j][i];
     54            if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.U8[j][i] & maskVal)) {
     55                smooth_im->data.F32[j][i] = 0.0;
     56            } else {
     57                smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i];
     58            }
     59        }
    6060    }
    6161    psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot"));
    6262
    63     // optionally save example images under trace 
     63    // optionally save example images under trace
    6464    if (psTraceGetLevel("psphot") > 5) {
    65         char name[64];
    66         sprintf (name, "snsmooth.v%d.fits", pass);
    67         psphotSaveImage (NULL, smooth_im, name);
     65        char name[64];
     66        sprintf (name, "snsmooth.v%d.fits", pass);
     67        psphotSaveImage (NULL, smooth_im, name);
    6868    }
    6969
    … …  
    7373    // signal/noise limit for the detected peaks
    7474    if (pass == 1) {
    75         NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
     75        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
    7676    } else {
    77         NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
    78     }   
    79    
     77        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
     78    }
     79
    8080    // we need to define the threshold based on the value of NSIGMA_PEAK and the applied smoothing
    8181    // gaussian SIGMA.  a peak in the significance image has an effective S/N for faint sources
    … …  
    8787    // terms of smooth_im peak counts Io, for a desired S/N limit corresponds to
    8888    // S/N = sqrt(Io)*4*pi*sigma_sm^2
    89     // thus, the threshold is: 
     89    // thus, the threshold is:
    9090    float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH);
    9191    if (effArea < 1) {
    92         effArea = 1;                    // never less than a pixel
     92        effArea = 1;                    // never less than a pixel
    9393    }
    9494    float threshold = PS_SQR(NSIGMA_PEAK) / effArea;
    … …  
    9797    psArray *peaks = pmFindImagePeaks (smooth_im, threshold);
    9898    if (peaks == NULL) {
    99         // XXX this may also be due to a programming or config error
    100         // XXX do we need to set something in the readout->analysis to indicate that
    101         // we tried and failed to find peaks (something in the header data)
    102         psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
    103         return false;
     99        // XXX this may also be due to a programming or config error
     100        // XXX do we need to set something in the readout->analysis to indicate that
     101        // we tried and failed to find peaks (something in the header data)
     102        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
     103        return false;
    104104    }
    105    
     105
    106106    // correct the peak values to S/N = sqrt(value*effArea)
    107107    // get the peak flux from the unsmoothed image
    … …  
    110110    int col0 = readout->image->col0;
    111111    for (int i = 0; i < peaks->n; i++) {
    112         pmPeak *peak = peaks->data[i];
    113         peak->SN = sqrt(peak->value*effArea);
    114         peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
     112        pmPeak *peak = peaks->data[i];
     113        peak->SN = sqrt(peak->value*effArea);
     114        peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
    115115    }
    116116
    117117    // limit the total number of returned peak as specified
    118118    if (pass == 1) {
    119         psArraySort (peaks, pmPeakSortBySN);
    120         int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX");
    121         if (NMAX && (peaks->n > NMAX)) {
    122             psArray *tmpPeaks = psArrayAllocEmpty (NMAX);
    123             for (int i = 0; i < NMAX; i++) {
    124                 psArrayAdd (tmpPeaks, 100, peaks->data[i]);
    125             }
    126             psFree (peaks);
    127             peaks = tmpPeaks;
    128         }
    129     }   
     119        psArraySort (peaks, pmPeakSortBySN);
     120        int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX");
     121        if (NMAX && (peaks->n > NMAX)) {
     122            psArray *tmpPeaks = psArrayAllocEmpty (NMAX);
     123            for (int i = 0; i < NMAX; i++) {
     124                psArrayAdd (tmpPeaks, 100, peaks->data[i]);
     125            }
     126            psFree (peaks);
     127            peaks = tmpPeaks;
     128        }
     129    }
    130130
    131131    // optional dump of all peak data
    132132    char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE");
    133133    if (status && (output != NULL) && (output[0])) {
    134         pmPeaksWriteText (peaks, output);
     134        pmPeaksWriteText (peaks, output);
    135135    }
    136136    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot"));
    … …  
    142142    // want to do that
    143143    //
    144     if (returnFootprints) {     // We want an array of pmFootprint, not pmPeak
    145         int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN");
    146         if (!status) {
    147             npixMin = 1;
    148         }
    149         float FOOTPRINT_NSIGMA_LIMIT =
    150             psMetadataLookupS32(&status, recipe,
    151                                 (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2");
    152         if (!status) {
    153             FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK;
    154         }
    155         threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea;
     144    if (returnFootprints) {     // We want an array of pmFootprint, not pmPeak
     145        int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN");
     146        if (!status) {
     147            npixMin = 1;
     148        }
     149        float FOOTPRINT_NSIGMA_LIMIT =
     150            psMetadataLookupS32(&status, recipe,
     151                                (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2");
     152        if (!status) {
     153            FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK;
     154        }
     155        threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea;
    156156
    157         psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin);
    158         pmPeaksAssignToFootprints(footprints, peaks);
     157        psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin);
     158        pmPeaksAssignToFootprints(footprints, peaks);
    159159
    160         psFree(peaks);
    161         peaks = footprints;             // well, you know what I mean
     160        psFree(peaks);
     161        peaks = footprints;             // well, you know what I mean
    162162    }
    163163
    … …  
    174174 */
    175175psErrorCode
    176 psphotCullPeaks(const psImage *image,   // the image wherein lives the footprint
    177                 const psImage *weight,  // corresponding variance image
    178                 const psMetadata *recipe,
    179                 psArray *footprints) {  // array of pmFootprints
     176psphotCullPeaks(const psImage *image,   // the image wherein lives the footprint
     177                const psImage *weight,  // corresponding variance image
     178                const psMetadata *recipe,
     179                psArray *footprints) {  // array of pmFootprints
    180180    bool status = false;
    181181    float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA");
    182182    if (!status) {
    183         nsigma_delta = 0; // min.
     183        nsigma_delta = 0; // min.
    184184    }
    185185    float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN");
    186186    if (!status) {
    187         nsigma_min = 0;
     187        nsigma_min = 0;
    188188    }
    189189    const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV");
    190190
    191191    return pmFootprintArrayCullPeaks(image, weight, footprints,
    192                                      nsigma_delta, nsigma_min*skyStdev);
     192                                     nsigma_delta, nsigma_min*skyStdev);
    193193}

  • trunk/psphot/src/psphotFitSet.c

    r13035 r13900  
    22
    33// This is not used in main psphot code (only in psphotModelTest.c)
    4 bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode) {
     4bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode, psMaskType maskVal) {
    55
    66    double x, y, Io;
    … …  
    2525
    2626    // XXX pmSourceFitSet must cache the modelFlux?
    27     pmSourceFitSet (source, modelSet, mode);
     27    pmSourceFitSet (source, modelSet, mode, maskVal);
    2828
    2929    // write out positive object
    … …  
    3333    for (int i = 0; i < modelSet->n; i++) {
    3434        pmModel *model = modelSet->data[i];
    35         pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     35        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    3636
    3737        fprintf (stderr, "output parameters (obj %d):\n", i);

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r13375 r13900  
    55// model with selected pixels, and the fit radius must be defined
    66
    7 // given the set of sources, each of which points to the pixels in the 
    8 // science image, we construct a set of simulated sources with their own pixels. 
    9 // these are used to determine the simultaneous linear fit of fluxes. 
     7// given the set of sources, each of which points to the pixels in the
     8// science image, we construct a set of simulated sources with their own pixels.
     9// these are used to determine the simultaneous linear fit of fluxes.
    1010// the analysis is performed wrt the simulated pixel values
    1111
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER);
    1313
    14 bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
     14bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final, psMaskType maskVal) {
    1515
    1616    bool status;
    … …  
    4141    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    4242    if (!status) {
    43         SKY_FIT_ORDER = 0;
     43        SKY_FIT_ORDER = 0;
    4444    }
    4545    bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
    4646    if (!status) {
    47         SKY_FIT_LINEAR = false;
     47        SKY_FIT_LINEAR = false;
    4848    }
    4949
    … …  
    5555        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    5656        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    57         if (source->type == PM_SOURCE_TYPE_STAR &&
     57        if (source->type == PM_SOURCE_TYPE_STAR &&
    5858            source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    5959        if (final) {
    … …  
    6767        y = source->peak->yf;
    6868
    69         // is the source in the region of interest?
     69        // is the source in the region of interest?
    7070        if (x < AnalysisRegion.x0) continue;
    7171        if (y < AnalysisRegion.y0) continue;
    … …  
    9797        psSparseMatrixElement (sparse, i, i, f);
    9898
    99         // the formal error depends on the weighting scheme
    100         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    101             float var = pmSourceModelDotModel (SRCi, SRCi, false);
    102             errors->data.F32[i] = 1.0 / sqrt(var);
    103         } else {
    104             errors->data.F32[i] = 1.0 / sqrt(f);
    105         }
     99        // the formal error depends on the weighting scheme
     100        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
     101            float var = pmSourceModelDotModel (SRCi, SRCi, false);
     102            errors->data.F32[i] = 1.0 / sqrt(var);
     103        } else {
     104            errors->data.F32[i] = 1.0 / sqrt(f);
     105        }
    106106
    107107
    … …  
    110110        psSparseVectorElement (sparse, i, f);
    111111
    112         // add the per-source weights (border region)
    113         switch (SKY_FIT_ORDER) {
    114           case 1:
    115             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
    116             psSparseBorderElementB (border, i, 1, f);
    117             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
    118             psSparseBorderElementB (border, i, 2, f);
    119 
    120           case 0:
    121             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
    122             psSparseBorderElementB (border, i, 0, f);
    123             break;
    124 
    125           default:
    126             psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER);
    127             break;
    128         }
     112        // add the per-source weights (border region)
     113        switch (SKY_FIT_ORDER) {
     114          case 1:
     115            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
     116            psSparseBorderElementB (border, i, 1, f);
     117            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
     118            psSparseBorderElementB (border, i, 2, f);
     119
     120          case 0:
     121            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
     122            psSparseBorderElementB (border, i, 0, f);
     123            break;
     124
     125          default:
     126            psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER);
     127            break;
     128        }
    129129
    130130        // loop over all other stars following this one
    … …  
    159159    psVector *skyfit = NULL;
    160160    if (SKY_FIT_LINEAR) {
    161         psSparseBorderSolve (&norm, &skyfit, constraint, border, 5);
    162         fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
     161        psSparseBorderSolve (&norm, &skyfit, constraint, border, 5);
     162        fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
    163163    } else {
    164         norm = psSparseSolve (NULL, constraint, sparse, 5);
    165         skyfit = NULL;
     164        norm = psSparseSolve (NULL, constraint, sparse, 5);
     165        skyfit = NULL;
    166166    }
    167167
    … …  
    169169    for (int i = 0; i < fitSources->n; i++) {
    170170        pmSource *source = fitSources->data[i];
    171         pmModel *model = pmSourceGetModel (NULL, source);
     171        pmModel *model = pmSourceGetModel (NULL, source);
    172172
    173173        // assign linearly-fitted normalization
    … …  
    177177        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
    178178        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
    179         // XXX is the value of 'errors' modified by the sky fit?
     179        // XXX is the value of 'errors' modified by the sky fit?
    180180
    181181        // subtract object
    182         pmSourceSub (source, PM_MODEL_OP_FULL);
     182        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    183183        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
    184184        if (!final) source->mode |= PM_SOURCE_MODE_TEMPSUB;
    185         // XXX not sure about the use of TEMPSUB
     185        // XXX not sure about the use of TEMPSUB
    186186    }
    187187
    … …  
    190190        pmSource *source = fitSources->data[i];
    191191        pmModel *model = pmSourceGetModel (NULL, source);
    192         pmSourceChisq (model, source->pixels, source->maskObj, source->weight);
     192        pmSourceChisq (model, source->pixels, source->maskObj, source->weight, maskVal);
    193193    }
    194194
    … …  
    219219    for (int i = 0; i < sources->n; i++) {
    220220        pmSource *source = sources->data[i];
    221         pmModel *model = pmSourceGetModel (NULL, source);
    222         if (model == NULL) continue;
     221        pmModel *model = pmSourceGetModel (NULL, source);
     222        if (model == NULL) continue;
    223223        float x = model->params->data.F32[PM_PAR_XPOS];
    224224        float y = model->params->data.F32[PM_PAR_YPOS];
    225         psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    226     }   
     225        psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
     226    }
    227227
    228228    // accumulate the image statistics from the masked regions
    … …  
    233233    double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
    234234    w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
    235    
     235
    236236    int col0 = readout->image->col0;
    237237    int row0 = readout->image->row0;
    238238
    239239    for (int j = 0; j < readout->image->numRows; j++) {
    240         for (int i = 0; i < readout->image->numCols; i++) {
    241             if (mask[j][i]) continue;
    242             if (constant_weights) {
    243                 wt = 1.0;
    244             } else {
    245                 wt = weight[j][i];
    246             }
    247             f = image[j][i];
    248             w   += 1/wt;
    249             fo  += f/wt;
    250             if (SKY_FIT_ORDER == 0) continue;
    251 
    252             xc  = i + col0;
    253             yc  = j + row0;
    254             x  +=    xc/wt;
    255             y  +=    yc/wt;
    256             x2 += xc*xc/wt;
    257             xy += xc*yc/wt;
    258             y2 += yc*yc/wt;
    259             fx +=  f*xc/wt;
    260             fy +=  f*yc/wt;
    261         }
     240        for (int i = 0; i < readout->image->numCols; i++) {
     241            if (mask[j][i]) continue;
     242            if (constant_weights) {
     243                wt = 1.0;
     244            } else {
     245                wt = weight[j][i];
     246            }
     247            f = image[j][i];
     248            w   += 1/wt;
     249            fo  += f/wt;
     250            if (SKY_FIT_ORDER == 0) continue;
     251
     252            xc  = i + col0;
     253            yc  = j + row0;
     254            x  +=    xc/wt;
     255            y  +=    yc/wt;
     256            x2 += xc*xc/wt;
     257            xy += xc*yc/wt;
     258            y2 += yc*yc/wt;
     259            fx +=  f*xc/wt;
     260            fy +=  f*yc/wt;
     261        }
    262262    }
    263263
    … …  
    269269    psSparseBorderElementT (border, 0, 0, w);
    270270    if (SKY_FIT_ORDER > 0) {
    271         psSparseBorderElementG (border, 0, fx);
    272         psSparseBorderElementG (border, 0, fy);
    273         psSparseBorderElementT (border, 1, 0, x);
    274         psSparseBorderElementT (border, 2, 0, y);
    275         psSparseBorderElementT (border, 0, 1, x);
    276         psSparseBorderElementT (border, 1, 1, x2);
    277         psSparseBorderElementT (border, 2, 1, xy);
    278         psSparseBorderElementT (border, 0, 2, y);
    279         psSparseBorderElementT (border, 1, 2, xy);
    280         psSparseBorderElementT (border, 2, 2, y2);
    281     }   
     271        psSparseBorderElementG (border, 0, fx);
     272        psSparseBorderElementG (border, 0, fy);
     273        psSparseBorderElementT (border, 1, 0, x);
     274        psSparseBorderElementT (border, 2, 0, y);
     275        psSparseBorderElementT (border, 0, 1, x);
     276        psSparseBorderElementT (border, 1, 1, x2);
     277        psSparseBorderElementT (border, 2, 1, xy);
     278        psSparseBorderElementT (border, 0, 2, y);
     279        psSparseBorderElementT (border, 1, 2, xy);
     280        psSparseBorderElementT (border, 2, 2, y2);
     281    }
    282282    return true;
    283283}

  • trunk/psphot/src/psphotGrowthCurve.c

    r13035 r13900  
    77// XXX add a option to turn off the curve-of-growth (ie, make the apMag = fitMag everywhere);
    88
    9 bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore) {
     9bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psMaskType maskVal) {
    1010
    1111    // bool status;
    … …  
    4545
    4646    // loop over a range of source fluxes
    47     // no need to interpolate since we have force the object center 
     47    // no need to interpolate since we have force the object center
    4848    // to 0.5, 0.5 above
    4949    for (int i = 0; i < psf->growth->radius->n; i++) {
    … …  
    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
     55        // NOTE: we use pmModelAdd not pmSourceAdd because we are not working with a normal source
    5656        psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK);
    57         pmModelAdd (image, mask, model, PM_MODEL_OP_FULL);
    58         pmSourcePhotometryAper (&apMag, model, image, mask);
     57        pmModelAdd (image, mask, model, PM_MODEL_OP_FULL, maskVal);
     58        pmSourcePhotometryAper (&apMag, model, image, mask, maskVal);
    5959        psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK));
    6060
    61         // the 'ignore' mode is for testing
    62         if (ignore) {
    63             psf->growth->apMag->data.F32[i] = fitMag;
    64         } else {
    65             psf->growth->apMag->data.F32[i] = apMag;
    66         }
     61        // the 'ignore' mode is for testing
     62        if (ignore) {
     63            psf->growth->apMag->data.F32[i] = fitMag;
     64        } else {
     65            psf->growth->apMag->data.F32[i] = apMag;
     66        }
    6767    }
    6868

  • trunk/psphot/src/psphotGuessModels.c

    r13569 r13900  
    1717
    1818// construct an initial PSF model for each object
    19 bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
     19bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) {
    2020
    2121  psTimerStart ("psphot");
    … …  
    5757    pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
    5858    if (modelPSF == NULL) {
    59         psError(PSPHOT_ERR_PSF, false,
    60                 "Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
    61                 source->peak->y, source->peak->x);
    62         //
    63         // Try the centre of the image
    64         //
    65         modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
    66         modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
    67         modelPSF = pmModelFromPSF (modelEXT, psf);
    68         if (modelPSF == NULL) {
    69             psError(PSPHOT_ERR_PSF, false,
    70                     "Failed to determine PSF model at centre of image");
    71             psFree(modelEXT);
    72             return false;
    73         }
     59        psError(PSPHOT_ERR_PSF, false,
     60                "Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
     61                source->peak->y, source->peak->x);
     62        //
     63        // Try the centre of the image
     64        //
     65        modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
     66        modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
     67        modelPSF = pmModelFromPSF (modelEXT, psf);
     68        if (modelPSF == NULL) {
     69            psError(PSPHOT_ERR_PSF, false,
     70                    "Failed to determine PSF model at centre of image");
     71            psFree(modelEXT);
     72            return false;
     73        }
    7474
    75         source->mode |= PM_SOURCE_MODE_BADPSF;
     75        source->mode |= PM_SOURCE_MODE_BADPSF;
    7676    }
    7777    psFree (modelEXT);
    … …  
    8585    source->modelPSF->residuals = psf->residuals;
    8686
    87     pmSourceCacheModel (source);
     87    pmSourceCacheModel (source, maskVal);
    8888  }
    8989  psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));

  • trunk/psphot/src/psphotImageMedian.c

    r13866 r13900  
    4040// generate the median in NxN boxes, clipping heavily
    4141// linear interpolation to generate full-scale model
    42 bool psphotImageMedian (pmConfig *config, const pmFPAview *view)
     42bool psphotImageMedian (pmConfig *config, const pmFPAview *view, psMaskType maskVal)
    4343{
    4444    bool status = true;
    … …  
    195195            // XXX don't bother trying if there are no valid pixels...
    196196
    197             if (psImageBackground(stats, subset, submask, 0xff, rng)) {
     197            if (psImageBackground(stats, subset, submask, maskVal, rng)) {
    198198                if (stats->options & PS_STAT_ROBUST_QUARTILE) {
    199199                    modelData[iy][ix] = stats->robustMedian;
    … …  
    205205                psStatsOptions currentOptions = stats->options;
    206206                stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
    207                 if (!psImageBackground(stats, subset, submask, 0xff, rng)) {
     207                if (!psImageBackground(stats, subset, submask, maskVal, rng)) {
    208208                    psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
    209209                               "(%dx%d, (row0,col0) = (%d,%d)",

  • trunk/psphot/src/psphotMagnitudes.c

    r12792 r13900  
    22
    33bool psphotMagnitudes(psArray *sources,
    4                       psMetadata *recipe,
    5                       pmPSF *psf,
    6                       pmReadout *background)
     4                      psMetadata *recipe,
     5                      pmPSF *psf,
     6                      pmReadout *background,
     7                      psMaskType maskVal,
     8                      psMaskType mark)
    79{
    810    bool status = false;
    … …  
    1315    pmSourceMagnitudesInit (recipe);
    1416
    15     // XXX require that we have a background model, or 
     17    // XXX require that we have a background model, or
    1618    // allow it to be missing, setting local sky to 0.0?
    1719    PS_ASSERT (background, false);
    1820
    19     // the binning details are saved on the analysis metadata 
     21    // the binning details are saved on the analysis metadata
    2022    psImageBinning *binning = psMetadataLookupPtr(&status, recipe, "PSPHOT.BACKGROUND.BINNING");
    2123    PS_ASSERT (status, false);
    … …  
    2931
    3032    for (int i = 0; i < sources->n; i++) {
    31         pmSource *source = (pmSource *) sources->data[i];
    32         status = pmSourceMagnitudes (source, psf, photMode);
    33         if (status) Nap ++;
     33        pmSource *source = (pmSource *) sources->data[i];
     34        status = pmSourceMagnitudes (source, psf, photMode, maskVal, mark);
     35        if (status) Nap ++;
    3436
    35         source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, background->image, binning);
    36         if (isnan(source->sky) && false) {
    37           psError(PSPHOT_ERR_SKY, false, "Setting pmSource.sky");
    38           psErrorStackPrint(NULL, " ");
    39           psErrorClear();
    40         }
    41     }   
     37        source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, background->image, binning);
     38        if (isnan(source->sky) && false) {
     39          psError(PSPHOT_ERR_SKY, false, "Setting pmSource.sky");
     40          psErrorStackPrint(NULL, " ");
     41          psErrorClear();
     42        }
     43    }
    4244
    4345    psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot"), sources->n, Nap);

  • trunk/psphot/src/psphotMakeResiduals.c

    r13806 r13900  
    11# include "psphotInternal.h"
    22
    3 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) {
     3bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) {
    44
    55    bool status, isPSF;
    … …  
    9494        psImage *weight = psImageCopy (NULL, source->weight,   PS_TYPE_F32);
    9595        psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_U8);
    96         pmModelSub (image, mask, model, PM_MODEL_OP_FUNC);
     96        pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal);
    9797
    9898        // re-normalize image and weight

  • trunk/psphot/src/psphotMaskReadout.c

    r13593 r13900  
    11# include "psphotInternal.h"
    22
    3 bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, pmConfig *config) {
     3bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, psMaskType maskVal) {
    44
    55    bool status;
    … …  
    1919
    2020    // psImageKeepRegion assumes the region refers to the parent coordinates
    21     psImageKeepRegion (readout->mask, keep, "OR", pmConfigMask("BAD", config));
     21    psImageKeepRegion (readout->mask, keep, "OR", maskVal);
    2222
    2323    return true;

  • trunk/psphot/src/psphotModelTest.c

    r13035 r13900  
    33
    44// XXX consider this function : add more test information?
    5 bool psphotModelTest (pmReadout *readout, psMetadata *recipe) {
     5bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {
    66
    77    bool status;
    … …  
    9494
    9595    // find the local sky
    96     status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     96    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark);
    9797    if (!status) psAbort("pmSourceLocalSky error");
    9898
    … …  
    159159
    160160    // define the pixels used for the fit
    161     psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", PM_MASK_MARK);
     161    psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", mark);
    162162
    163163    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    164164    if (status) {
    165         status = psphotFitSet (source, model, fitset, fitMode);
     165        status = psphotFitSet (source, model, fitset, fitMode, maskVal);
    166166        exit (0);
    167167    }
    168168
    169     status = pmSourceFitModel (source, model, fitMode);
     169    status = pmSourceFitModel (source, model, fitMode, maskVal);
    170170
    171171    // measure the source mags
    172172    pmSourcePhotometryModel (&fitMag, model);
    173     pmSourcePhotometryAper  (&obsMag, model, source->pixels, source->maskObj);
     173    pmSourcePhotometryAper  (&obsMag, model, source->pixels, source->maskObj, maskVal);
    174174    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
    175175
    … …  
    178178
    179179    // subtract object, leave local sky
    180     pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     180    pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    181181
    182182    fprintf (stderr, "output parameters: \n");

  • trunk/psphot/src/psphotMosaicChip.c

    r12792 r13900  
    2828    psTrace("pmChipMosaic", 5, "mosaic chip %s to %s (xbin,ybin: %d,%d to %d,%d)\n",
    2929            in->name, out->name, in->xBin, in->yBin, out->xBin, out->yBin);
    30     status = pmChipMosaic(outChip, inChip, true);
     30    status = pmChipMosaic(outChip, inChip, true, pmConfigMask("BLANK", config));
    3131    return status;
    3232}

  • trunk/psphot/src/psphotReadout.c

    r13804 r13900  
    1212    }
    1313
     14    // Interpret the mask values
     15    const char *maskValStr = psMetadataLookupStr(NULL, recipe, "MASKVAL"); // String with mask names
     16    if (!maskValStr) {
     17        psError(PSPHOT_ERR_CONFIG, false, "Missing recipe item: MASKVAL(STR)");
     18        return false;
     19    }
     20    psMaskType maskVal = pmConfigMask(maskValStr, config); // Mask values to mask against
     21    psMaskType maskMark = pmConfigMask("MARK", config); // Mask value for marking
     22    psMaskType maskSat = pmConfigMask("SAT", config); // Mask value for saturated pixels
     23    psMaskType maskBad = pmConfigMask("BAD", config); // Mask value for bad pixels
     24
    1425    // find the currently selected readout
    1526    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    … …  
    2435
    2536    // generate mask & weight images if they don't already exit
    26     psMaskType satMask = pmConfigMask("SAT", config);
    27     psMaskType badMask = pmConfigMask("BAD", config);
    2837    if (!readout->mask) {
    29         if (!pmReadoutGenerateMask(readout, satMask, badMask)) {
    30             psError (PSPHOT_ERR_CONFIG, false, "trouble creating mask");
    31             return false;
    32         }
     38        if (!pmReadoutGenerateMask(readout, maskSat, maskBad)) {
     39            psError (PSPHOT_ERR_CONFIG, false, "trouble creating mask");
     40            return false;
     41        }
    3342    }
    3443    if (!readout->weight) {
    35         if (!pmReadoutGenerateWeight(readout, true)) {
    36             psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight");
    37             return false;
    38         }
     44        if (!pmReadoutGenerateWeight(readout, true)) {
     45            psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight");
     46            return false;
     47        }
    3948    }
    4049
    … …  
    4655
    4756    // I have a valid mask, now mask in the analysis region of interest
    48     psphotMaskReadout (readout, recipe, config);
     57    psphotMaskReadout (readout, recipe, maskBad);
    4958
    5059    // run a single-model test if desired
    51     psphotModelTest (readout, recipe);
     60    psphotModelTest (readout, recipe, maskVal, maskMark);
    5261
    5362    if (psTraceGetLevel("psphot") >= 5) {
    … …  
    6271
    6372    // generate a background model (median, smoothed image)
    64     if (!psphotImageMedian (config, view)) {
     73    if (!psphotImageMedian (config, view, maskVal)) {
    6574        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL);
    6675    }
    … …  
    7483    psArray *footprints = NULL;
    7584    if (useFootprints) {
    76        footprints = psphotFindPeaks (readout, recipe, useFootprints, 1);
     85       footprints = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal);
    7786       int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS");
    7887       if (growRadius > 0) {
    … …  
    8594       peaks = pmFootprintArrayToPeaks(footprints);
    8695    } else {
    87        peaks = psphotFindPeaks (readout, recipe, useFootprints, 1);
     96       peaks = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal);
    8897    }
    8998
    … …  
    97106
    98107    // construct sources and measure basic stats
    99     psArray *sources = psphotSourceStats (readout, recipe, peaks);
     108    psArray *sources = psphotSourceStats (readout, recipe, peaks, maskVal, maskMark);
    100109    if (!sources) return false;
    101110    psFree (peaks);
    … …  
    114123
    115124    // classify sources based on moments, brightness
    116     if (!psphotRoughClass (sources, recipe, true)) {
     125    if (!psphotRoughClass (sources, recipe, true, maskSat)) {
    117126        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    118127        return psphotReadoutCleanup (config, readout, recipe, NULL, sources);
    … …  
    126135    if (!psf) {
    127136        // use bright stellar objects to measure PSF
    128         psf = psphotChoosePSF (readout, sources, recipe);
     137        psf = psphotChoosePSF (readout, sources, recipe, maskVal, maskMark);
    129138        if (psf == NULL) {
    130139            psLogMsg ("psphot", 3, "failure to construct a psf model");
    … …  
    143152
    144153    // construct an initial model for each object
    145     psphotGuessModels (readout, sources, recipe, psf);
     154    psphotGuessModels (readout, sources, recipe, psf, maskVal);
    146155
    147156    // XXX test output of models
    … …  
    151160
    152161    // linear PSF fit to peaks
    153     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
     162    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE, maskVal);
    154163    if (dump) psphotSaveImage (NULL, readout->image,  "image.v1.fits");
    155164
    … …  
    162171
    163172    // non-linear PSF and EXT fit to brighter sources
    164     psphotBlendFit (readout, sources, recipe, psf);
     173    psphotBlendFit (readout, sources, recipe, psf, maskVal);
    165174    if (dump) psphotSaveImage (NULL, readout->image,  "image.v2.fits");
    166175
    167176    // replace all sources
    168     psphotReplaceAll (sources);
     177    psphotReplaceAll (sources, maskVal);
    169178    if (dump) psphotSaveImage (NULL, readout->image,  "image.v3.fits");
    170179
    171180    // linear PSF fit to remaining peaks
    172     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
     181    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal);
    173182    if (dump) psphotSaveImage (NULL, readout->image,  "image.v4.fits");
    174183    if (!strcasecmp (breakPt, "PASS1")) {
    … …  
    183192
    184193    // add noise for subtracted objects
    185     psphotAddNoise (readout, sources, recipe, true);
     194    psphotAddNoise (readout, sources, recipe, true, maskVal);
    186195
    187196    // find the peaks in the image
    188197    psArray *newPeaks;
    189198    if (useFootprints) {
    190        psArray *newFootprints = psphotFindPeaks (readout, recipe, useFootprints, 2);
     199       psArray *newFootprints = psphotFindPeaks (readout, recipe, useFootprints, 2, maskVal);
    191200
    192201       int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS_2");
    … …  
    209218       psFree(mergedFootprints);
    210219    } else {
    211        newPeaks = psphotFindPeaks(readout, recipe, useFootprints, 2);
     220       newPeaks = psphotFindPeaks(readout, recipe, useFootprints, 2, maskVal);
    212221    }
    213222
    214223    // remove noise for subtracted objects
    215     psphotAddNoise (readout, sources, recipe, false);
     224    psphotAddNoise (readout, sources, recipe, false, maskVal);
    216225
    217226    // define new sources based on the new peaks
    218     psArray *newSources = psphotSourceStats (readout, recipe, newPeaks);
     227    psArray *newSources = psphotSourceStats (readout, recipe, newPeaks, maskVal, maskMark);
    219228    psFree (newPeaks);
    220229
    221230    // set source type
    222     psphotRoughClass (newSources, recipe, false);
     231    psphotRoughClass (newSources, recipe, false, maskSat);
    223232
    224233    // create full input models
    225     psphotGuessModels (readout, newSources, recipe, psf);
     234    psphotGuessModels (readout, newSources, recipe, psf, maskVal);
    226235
    227236    // replace all sources
    228     psphotReplaceAll (sources);
     237    psphotReplaceAll (sources, maskVal);
    229238    if (dump) psphotSaveImage (NULL, readout->image,  "image.v5.fits");
    230239
    … …  
    234243
    235244    // linear PSF fit to remaining peaks
    236     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
     245    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal);
    237246    if (dump) psphotSaveImage (NULL, readout->image,  "image.v6.fits");
    238247
    … …  
    243252
    244253    // measure aperture photometry corrections
    245     if (!psphotApResid (readout, sources, recipe, psf)) {
     254    if (!psphotApResid (readout, sources, recipe, psf, maskVal, maskMark)) {
    246255        psTrace ("psphot", 4, "failure on psphotApResid");
    247256        psError(PSPHOT_ERR_PHOTOM, false, "Measure aperture photometry corrections");
    … …  
    251260    // calculate source magnitudes
    252261    pmReadout *background = psphotSelectBackground (config, view, false);
    253     psphotMagnitudes(sources, recipe, psf, background);
     262    psphotMagnitudes(sources, recipe, psf, background, maskVal, maskMark);
    254263
    255264    // replace failed sources?

  • trunk/psphot/src/psphotReplaceUnfit.c

    r13035 r13900  
    22
    33// replace the flux for sources which failed
    4 bool psphotReplaceUnfit (psArray *sources) {
     4bool psphotReplaceUnfit (psArray *sources, psMaskType maskVal) {
    55
    66    pmSource *source;
    … …  
    1616
    1717    replace:
    18         pmSourceAdd (source, PM_MODEL_OP_FULL);
    19         source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
    20         source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
     18        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     19        source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
     20        source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
    2121    }
    2222    psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%ld objects)\n", psTimerMark ("psphot"), sources->n);
    … …  
    2424}
    2525
    26 bool psphotReplaceAll (psArray *sources) {
     26bool psphotReplaceAll (psArray *sources, psMaskType maskVal) {
    2727
    2828    pmSource *source;
    … …  
    3636      if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
    3737
    38       pmSourceAdd (source, PM_MODEL_OP_FULL);
     38      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    3939      source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
    4040    }
    … …  
    4444
    4545// add source, if the source has been subtracted (or if we ignore the state)
    46 bool psphotAddWithTest (pmSource *source, bool useState) {
     46bool psphotAddWithTest (pmSource *source, bool useState, psMaskType maskVal) {
    4747
    4848    // what is current state? (true : add; false : sub)
    4949    bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
    5050    if (state && useState) return true;
    51    
     51
    5252    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
    5353    if (!state || !useState) {
    54         pmSourceAdd (source, PM_MODEL_OP_FULL);
     54        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    5555    }
    5656    return true;
    … …  
    5858
    5959// sub source, if the source has been added (or if we ignore the state)
    60 bool psphotSubWithTest (pmSource *source, bool useState) {
     60bool psphotSubWithTest (pmSource *source, bool useState, psMaskType maskVal) {
    6161
    6262    // what is current state? (true : sub; false : add)
    … …  
    6666    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
    6767    if (!state || !useState) {
    68         pmSourceSub (source, PM_MODEL_OP_FULL);
     68        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    6969    }
    7070    return true;
    7171}
    7272
    73 // add or sub source replace or if the source has 
    74 bool psphotSetState (pmSource *source, bool curState) {
     73// add or sub source replace or if the source has
     74bool psphotSetState (pmSource *source, bool curState, psMaskType maskVal) {
    7575
    7676    // what is desired state? (true : add; false : sub)
    7777    bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
    7878    if (curState == newState) return true;
    79    
     79
    8080    if (curState && !newState) {
    81         pmSourceSub (source, PM_MODEL_OP_FULL);
     81        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    8282    }
    8383    if (newState && !curState) {
    84         pmSourceAdd (source, PM_MODEL_OP_FULL);
     84        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    8585    }
    8686    return true;

  • trunk/psphot/src/psphotRoughClass.c

    r13374 r13900  
    22
    33// 2006.02.02 : no leaks
    4 bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump) {
     4bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump, psMaskType maskSat) {
    55
    66    static pmPSFClump   psfClump;
    … …  
    1010
    1111    if (findPsfClump) {
    12         psfClump = pmSourcePSFClump (sources, recipe);
     12        psfClump = pmSourcePSFClump (sources, recipe);
    1313    } else if (!havePsfClump) {
    14         psError(PSPHOT_ERR_PROG, false, "You must find the PSF clump before reusing it");
     14        psError(PSPHOT_ERR_PROG, false, "You must find the PSF clump before reusing it");
    1515    } else {
    16         ;
     16        ;
    1717    }
    18    
    19     havePsfClump = false;               // we don't know if it's valid
     18
     19    havePsfClump = false;               // we don't know if it's valid
    2020    if (psfClump.X < 0) {
    21         psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump");
    22         return false;
     21        psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump");
     22        return false;
    2323    }
    2424    if (!psfClump.X || !psfClump.Y) {
    25         psError(PSPHOT_ERR_DATA, true, "Failed to find a valid PSF clump");
    26         return false;
     25        psError(PSPHOT_ERR_DATA, true, "Failed to find a valid PSF clump");
     26        return false;
    2727    }
    2828    psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
    … …  
    3030
    3131    // group into STAR, COSMIC, EXTENDED, SATURATED, etc.
    32     if (!pmSourceRoughClass (sources, recipe, psfClump)) {
    33         psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass");
    34         return false;
     32    if (!pmSourceRoughClass (sources, recipe, psfClump, maskSat)) {
     33        psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass");
     34        return false;
    3535    }
    3636
    … …  
    4040    psLogMsg ("psphot.roughclass", PS_LOG_INFO, "rough classification: %f sec\n", psTimerMark ("psphot"));
    4141
    42     havePsfClump = true;                        // we have set psfClump
    43    
     42    havePsfClump = true;                        // we have set psfClump
     43
    4444    return true;
    4545}

  • trunk/psphot/src/psphotSourceFits.c

    r13035 r13900  
    1616bool psphotFitSummary () {
    1717
    18     fprintf (stderr, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", 
    19              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
    20     return true;
    21 }
    22 
    23 bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) {
     18    fprintf (stderr, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
     19             NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
     20    return true;
     21}
     22
     23bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal) {
    2424
    2525    float x, y, dR;
    … …  
    2727    // if this source is not a possible blend, just fit as PSF
    2828    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
    29         bool status = psphotFitPSF (readout, source, psf);
     29        bool status = psphotFitPSF (readout, source, psf, maskVal);
    3030        return status;
    3131    }
    … …  
    6565        model->params->data.F32[PM_PAR_YPOS] = blend->peak->yf;
    6666
    67         // these should never be invalid values
    68         // XXX drop these tests eventually
     67        // these should never be invalid values
     68        // XXX drop these tests eventually
    6969        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");
    … …  
    8383
    8484    // fit PSF model (set/unset the pixel mask)
    85     pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
     85    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
    8686
    8787    // correct model chisq for flux trend
    … …  
    109109        psTrace ("psphot", 5, "fitted blend as PSF\n");
    110110
    111         // build cached model and subtract
    112         pmSourceCacheModel (blend);
    113         pmSourceSub (blend, PM_MODEL_OP_FULL);
     111        // build cached model and subtract
     112        pmSourceCacheModel (blend, maskVal);
     113        pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal);
    114114        blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    115115        blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
    … …  
    121121        psTrace ("psphot", 4, "failed on blend 0 of %ld\n", modelSet->n);
    122122        psFree (PSF);
    123         psFree (modelSet);
    124         psFree (sourceSet);
     123        psFree (modelSet);
     124        psFree (sourceSet);
    125125        return false;
    126126    }
    … …  
    134134
    135135    // build cached model and subtract
    136     pmSourceCacheModel (source);
    137     pmSourceSub (source, PM_MODEL_OP_FULL);
     136    pmSourceCacheModel (source, maskVal);
     137    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    138138    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    139139    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
    … …  
    141141}
    142142
    143 bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) {
     143bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal) {
    144144
    145145    double chiTrend;
    … …  
    155155
    156156    // fit PSF model (set/unset the pixel mask)
    157     pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF);
     157    pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal);
    158158
    159159    // correct model chisq for flux trend
    … …  
    173173
    174174    // build cached model and subtract
    175     pmSourceCacheModel (source);
    176     pmSourceSub (source, PM_MODEL_OP_FULL);
     175    pmSourceCacheModel (source, maskVal);
     176    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    177177    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
    178178    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
    … …  
    201201}
    202202
    203 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) {
     203bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal) {
    204204
    205205    bool okEXT, okDBL;
    … …  
    222222    pmSource *tmpSrc = pmSourceAlloc ();
    223223
    224     pmModel *EXT = psphotFitEXT (readout, source);
     224    pmModel *EXT = psphotFitEXT (readout, source, maskVal);
    225225    okEXT = psphotEvalEXT (tmpSrc, EXT);
    226226    chiEXT = EXT->chisq / EXT->nDOF;
    227227
    228     psArray *DBL = psphotFitDBL (readout, source);
     228    psArray *DBL = psphotFitDBL (readout, source, maskVal);
    229229    okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
    230230    okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
    … …  
    269269
    270270    // build cached model and subtract
    271     pmSourceCacheModel (source);
    272     pmSourceSub (source, PM_MODEL_OP_FULL);
     271    pmSourceCacheModel (source, maskVal);
     272    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    273273    psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
    274274
    … …  
    294294
    295295    // 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);
     296    pmSourceCacheModel (source, maskVal);
     297    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     298    pmSourceCacheModel (newSrc, maskVal);
     299    pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal);
    300300    psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
    301301
    … …  
    307307
    308308// fit a double PSF source to an extended blob
    309 psArray *psphotFitDBL (pmReadout *readout, pmSource *source) {
     309psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal) {
    310310
    311311    float dx, dy;
    … …  
    348348
    349349    // fit PSF model (set/unset the pixel mask)
    350     pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
     350    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
    351351    return (modelSet);
    352352}
    353353
    354 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) {
     354pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal) {
    355355
    356356    NfitEXT ++;
    … …  
    359359    pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT);
    360360    PS_ASSERT (EXT, NULL);
    361        
     361
    362362    // if (isnan(EXT->params->data.F32[1])) psAbort("nan in ext fit");
    363363
    … …  
    369369
    370370    // fit EXT (not PSF) model (set/unset the pixel mask)
    371     pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT);
     371    pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT, maskVal);
    372372    return (EXT);
    373373}

  • trunk/psphot/src/psphotSourcePlots.c

    r13006 r13900  
    11# include "psphotInternal.h"
    22
    3 bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe) {
     3bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal) {
    44
    55    // bool status = false;
    … …  
    1717
    1818    // counters to track the size of the image and area used in a row
    19     int dX = 0;                         // starting corner of next box
    20     int dY = 0;                         // height of row so far
    21     int NX = 20*DX;                     // full width of output image
    22     int NY = 0;                         // total height of output image
     19    int dX = 0;                         // starting corner of next box
     20    int dY = 0;                         // height of row so far
     21    int NX = 20*DX;                     // full width of output image
     22    int NY = 0;                         // total height of output image
    2323
    2424    // first, examine the PSF and SAT stars:
    … …  
    2828        pmSource *source = sources->data[i];
    2929
    30         bool keep = false;
     30        bool keep = false;
    3131        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
    3232        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
    33         if (!keep) continue;
     33        if (!keep) continue;
    3434
    35         // how does this subimage get placed into the output image?
    36         // DX = source->pixels->numCols
    37         // DY = source->pixels->numRows
     35        // how does this subimage get placed into the output image?
     36        // DX = source->pixels->numCols
     37        // DY = source->pixels->numRows
    3838
    39         if (dX + DX > NX) {
    40             // too wide for the rest of this row
    41             if (dX == 0) {
    42                 // alone on this row
    43                 NY += DY;
    44                 dX = 0;
    45                 dY = 0;
    46             } else {
    47                 // start the next row
    48                 NY += dY;
    49                 dX = DX;
    50                 dY = DY;
    51             }
    52         } else {
    53             // extend this row
    54             dX += DX;
    55             dY = PS_MAX (dY, DY);
    56         }
     39        if (dX + DX > NX) {
     40            // too wide for the rest of this row
     41            if (dX == 0) {
     42                // alone on this row
     43                NY += DY;
     44                dX = 0;
     45                dY = 0;
     46            } else {
     47                // start the next row
     48                NY += dY;
     49                dX = DX;
     50                dY = DY;
     51            }
     52        } else {
     53            // extend this row
     54            dX += DX;
     55            dY = PS_MAX (dY, DY);
     56        }
    5757    }
    5858
    … …  
    6161    psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32);
    6262
    63     int Xo = 0;                         // starting corner of next box
    64     int Yo = 0;                         // starting corner of next box
    65     dY = 0;                             // height of row so far
     63    int Xo = 0;                         // starting corner of next box
     64    int Yo = 0;                         // starting corner of next box
     65    dY = 0;                             // height of row so far
    6666
    6767    int nPSF = 0;
    6868    int nSAT = 0;
    69     int kapa = 0;                       // file descriptor for plotting routine
     69    int kapa = 0;                       // file descriptor for plotting routine
    7070
    7171    // first, examine the PSF and SAT stars:
    … …  
    7676        pmSource *source = sources->data[i];
    7777
    78         bool keep = false;
     78        bool keep = false;
    7979        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
    80             nPSF ++;
    81             keep = true;
    82         }
     80            nPSF ++;
     81            keep = true;
     82        }
    8383        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    84             nSAT ++;
    85             keep = true;
    86         }           
    87         if (!keep) continue;
     84            nSAT ++;
     85            keep = true;
     86        }
     87        if (!keep) continue;
    8888
    89         // how does this subimage get placed into the output image?
    90         // DX = source->pixels->numCols
    91         // DY = source->pixels->numRows
     89        // how does this subimage get placed into the output image?
     90        // DX = source->pixels->numCols
     91        // DY = source->pixels->numRows
    9292
    93         if (Xo + DX > NX) {
    94             // too wide for the rest of this row
    95             if (Xo == 0) {
    96                 // place source alone on this row
    97                 psphotAddWithTest (source, true); // replace source if subtracted
    98                 psphotRadialPlot (&kapa, "radial.plots.ps", source);
    99                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
     93        if (Xo + DX > NX) {
     94            // too wide for the rest of this row
     95            if (Xo == 0) {
     96                // place source alone on this row
     97                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     98                psphotRadialPlot (&kapa, "radial.plots.ps", source);
     99                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
    100100
    101                 psphotSubWithTest (source, false); // remove source (force)
    102                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
     101                psphotSubWithTest (source, false, maskVal); // remove source (force)
     102                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
    103103
    104                 psphotSetState (source, false); // replace source (has been subtracted)
    105                 Yo += DY;
    106                 Xo = 0;
    107                 dY = 0;
    108             } else {
    109                 // start the next row
    110                 Yo += dY;
    111                 Xo = 0;
    112                 psphotAddWithTest (source, true); // replace source if subtracted
    113                 psphotRadialPlot (&kapa, "radial.plots.ps", source);
    114                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
     104                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     105                Yo += DY;
     106                Xo = 0;
     107                dY = 0;
     108            } else {
     109                // start the next row
     110                Yo += dY;
     111                Xo = 0;
     112                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     113                psphotRadialPlot (&kapa, "radial.plots.ps", source);
     114                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
    115115
    116                 psphotSubWithTest (source, false); // remove source (force)
    117                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
    118                 psphotSetState (source, false); // replace source (has been subtracted)
     116                psphotSubWithTest (source, false, maskVal); // remove source (force)
     117                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
     118                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
    119119
    120                 Xo = DX;
    121                 dY = DY;
    122             }
    123         } else {
    124             // extend this row
    125             psphotAddWithTest (source, true); // replace source if subtracted
    126             psphotRadialPlot (&kapa, "radial.plots.ps", source);
    127             psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
     120                Xo = DX;
     121                dY = DY;
     122            }
     123        } else {
     124            // extend this row
     125            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     126            psphotRadialPlot (&kapa, "radial.plots.ps", source);
     127            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
    128128
    129             psphotSubWithTest (source, false); // remove source (force)
    130             psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
    131             psphotSetState (source, false); // replace source (has been subtracted)
     129            psphotSubWithTest (source, false, maskVal); // remove source (force)
     130            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
     131            psphotSetState (source, false, maskVal); // replace source (has been subtracted)
    132132
    133             Xo += DX;
    134             dY = PS_MAX (dY, DY);
    135         }
     133            Xo += DX;
     134            dY = PS_MAX (dY, DY);
     135        }
    136136    }
    137137

  • trunk/psphot/src/psphotSourceStats.c

    r12792 r13900  
    11# include "psphotInternal.h"
    22
    3 psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks)
     3psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks, psMaskType maskVal, psMaskType mark)
    44{
    55    bool     status  = false;
    … …  
    4040
    4141        // skip faint sources
    42         if (source->peak->SN < MIN_SN) {
     42        if (source->peak->SN < MIN_SN) {
    4343            psArrayAdd (sources, 100, source);
    4444            psFree (source);
    45             continue;
    46         }
     45            continue;
     46        }
    4747
    4848        // measure a local sky value
    4949        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
    50         status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     50        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark);
    5151        if (!status) {
    5252          psFree (source);
    53           Nfail ++;
     53          Nfail ++;
    5454          continue;
    5555        }
    … …  
    5757        // measure the local sky variance (needed if noise is not sqrt(signal))
    5858        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
    59         status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     59        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark);
    6060        if (!status) {
    6161          psFree (source);
    62           Nfail ++;
     62          Nfail ++;
    6363          continue;
    6464        }
    … …  
    7070            psArrayAdd (sources, 100, source);
    7171            psFree (source);
    72             Nmoments ++;
     72            Nmoments ++;
    7373            continue;
    7474        }
    … …  
    8383            psArrayAdd (sources, 100, source);
    8484            psFree (source);
    85             Nmoments ++;
     85            Nmoments ++;
    8686            continue;
    8787        }
    8888
    8989        psFree (source);
    90         Nfail ++;
     90        Nfail ++;
    9191        continue;
    9292    }
Note: See TracChangeset for help on using the changeset viewer.