IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25852


Ignore:
Timestamp:
Oct 15, 2009, 12:10:16 PM (17 years ago)
Author:
Paul Price
Message:

Putting non-writable PS_DATA_REGION type on the recipe was causing the configuration dumping to fail. pmReadout.analysis is a better place for these analysis data.

Location:
trunk/psphot/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphot.h

    r25755 r25852  
    182182bool psphotVisualShowFootprints (pmDetections *detections);
    183183bool psphotVisualShowMoments (psArray *sources);
    184 bool psphotVisualPlotMoments (psMetadata *recipe, psArray *sources);
     184bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    185185bool psphotVisualShowRoughClass (psArray *sources);
    186186bool psphotVisualShowPSFStars (psMetadata *recipe, pmPSF *psf, psArray *sources);
    … …  
    193193bool psphotVisualShowResidualImage (pmReadout *readout);
    194194bool psphotVisualPlotApResid (psArray *sources, float mean, float error);
    195 bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources);
     195bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    196196bool psphotVisualShowPetrosians (psArray *sources);
    197197
    … …  
    212212// bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef);
    213213// bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
    214 // bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, 
    215 //                               psVector *refRadius, psVector *meanSB,
    216 //                               psVector *petRatio, psVector *petRatioErr, psVector *fluxSum,
    217 //                               float petRadius, float ratioForRadius,
    218 //                               float petFlux, float radiusForFlux);
     214// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
     215//                               psVector *refRadius, psVector *meanSB,
     216//                               psVector *petRatio, psVector *petRatioErr, psVector *fluxSum,
     217//                               float petRadius, float ratioForRadius,
     218//                               float petFlux, float radiusForFlux);
    219219
    220220bool psphotImageQuality (psMetadata *recipe, psArray *sources);

  • trunk/psphot/src/psphotRoughClass.c

    r25755 r25852  
    11# include "psphotInternal.h"
    22
    3 # define CHECK_STATUS(S,MSG) {                                          \
    4         if (!status) {                                                  \
    5             psError(PSPHOT_ERR_CONFIG, false, "missing PSF Clump entry: %s\n", MSG); \
    6             return false;                                               \
    7         } }
     3# define CHECK_STATUS(S,MSG) {                                          \
     4        if (!status) {                                                  \
     5            psError(PSPHOT_ERR_CONFIG, false, "missing PSF Clump entry: %s\n", MSG); \
     6            return false;                                               \
     7        } }
    88
    9 bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *recipe, const bool havePSF);
     9bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF);
    1010
    1111// 2006.02.02 : no leaks
    … …  
    2424    int nRegion = 0;
    2525    for (int ix = 0; ix < NX; ix ++) {
    26         for (int iy = 0; iy < NY; iy ++) {
     26        for (int iy = 0; iy < NY; iy ++) {
    2727
    28             psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY);
    29             if (!psphotRoughClassRegion (nRegion, region, sources, recipe, havePSF)) {
    30                 psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n",
    31                         region->x0, region->y0, region->x1, region->y1);
    32                 psFree (region);
    33                 continue;
    34             }
    35             psFree (region);
    36             nRegion ++;
    37         }
     28            psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY);
     29            if (!psphotRoughClassRegion (nRegion, region, sources, readout->analysis, recipe, havePSF)) {
     30                psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n",
     31                        region->x0, region->y0, region->x1, region->y1);
     32                psFree (region);
     33                continue;
     34            }
     35            psFree (region);
     36            nRegion ++;
     37        }
    3838    }
    3939    psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegion);
    … …  
    4444    psLogMsg ("psphot.roughclass", PS_LOG_INFO, "rough classification: %f sec\n", psTimerMark ("psphot.rough"));
    4545
    46     psphotVisualPlotMoments (recipe, sources);
     46    psphotVisualPlotMoments (recipe, readout->analysis, sources);
    4747    psphotVisualShowRoughClass (sources);
    4848    // XXX better visualization: psphotVisualShowFlags (sources);
    … …  
    5151}
    5252
    53 bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *recipe, const bool havePSF) {
     53bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF) {
    5454
    5555    bool status;
    … …  
    6262
    6363    snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", nRegion);
    64     psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     64    psMetadata *regionMD = psMetadataLookupPtr (&status, target, regionName);
    6565    if (!regionMD) {
    66         // allocate the region metadata folder and add this region to it.
    67         regionMD = psMetadataAlloc();
    68         psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
    69         psFree (regionMD);
     66        // allocate the region metadata folder and add this region to it.
     67        regionMD = psMetadataAlloc();
     68        psMetadataAddMetadata (target, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
     69        psFree (regionMD);
    7070    }
    7171    psMetadataAddPtr (regionMD, PS_LIST_TAIL, "REGION", PS_DATA_REGION | PS_META_REPLACE, "psf clump region", region);
    7272
    7373    if (!havePSF) {
    74         // determine the PSF parameters from the source moment values
    75         // XXX why not save the psfClump as a PTR?
    76         psfClump = pmSourcePSFClump (region, sources, recipe);
    77         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
    78         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
    79         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
    80         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
     74        // determine the PSF parameters from the source moment values
     75        // XXX why not save the psfClump as a PTR?
     76        psfClump = pmSourcePSFClump (region, sources, target);
     77        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
     78        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
     79        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
     80        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
    8181    } else {
    82         // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
    83         psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    84         if (!status) {
    85             psLogMsg ("psphot", 4, "No PSF clump defined for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    86             return false;
    87         }           
    88         psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
    89         psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
    90         psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
     82        // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
     83        psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     84        if (!status) {
     85            psLogMsg ("psphot", 4, "No PSF clump defined for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     86            return false;
     87        }
     88        psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
     89        psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
     90        psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
    9191    }
    9292
    9393    if (psfClump.X < 0) {
    94         psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump");
    95         return false;
     94        psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump");
     95        return false;
    9696    }
    9797    if (!psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
    98         psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    99         return false;
     98        psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     99        return false;
    100100    }
    101101    psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
    … …  
    104104    // group into STAR, COSMIC, EXTENDED, SATURATED, etc.
    105105    if (!pmSourceRoughClass (region, sources, recipe, psfClump, maskSat)) {
    106         psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass");
    107         return false;
     106        psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass");
     107        return false;
    108108    }
    109109

  • trunk/psphot/src/psphotSourceSize.c

    r25755 r25852  
    8383    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size"));
    8484
    85     psphotVisualPlotSourceSize (recipe, sources);
     85    psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    8686    psphotVisualShowSourceSize (readout, sources);
    8787    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
    … …  
    103103    pmFootprint *footprint = peak->footprint;
    104104    if (!footprint) {
    105         // if we have not footprint, use the old code to mask by isophot
    106         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    107         return true;
     105        // if we have not footprint, use the old code to mask by isophot
     106        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     107        return true;
    108108    }
    109109
    110110    if (!footprint->spans) {
    111         // if we have no footprint, use the old code to mask by isophot
    112         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    113         return true;
     111        // if we have no footprint, use the old code to mask by isophot
     112        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     113        return true;
    114114    }
    115115
    116116    // mask all of the pixels covered by the spans of the footprint
    117117    for (int j = 1; j < footprint->spans->n; j++) {
    118         pmSpan *span1 = footprint->spans->data[j];
    119 
    120         int iy = span1->y;
    121         int xs = span1->x0;
    122         int xe = span1->x1;
    123 
    124         for (int ix = xs; ix < xe; ix++) {
    125             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    126         }
     118        pmSpan *span1 = footprint->spans->data[j];
     119
     120        int iy = span1->y;
     121        int xs = span1->x0;
     122        int xe = span1->x1;
     123
     124        for (int ix = xs; ix < xe; ix++) {
     125            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     126        }
    127127    }
    128128    return true;
    … …  
    147147    // mark the pixels in this row to the left, then the right
    148148    for (int ix = xo; ix >= 0; ix--) {
    149         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    150         if (SN > SN_LIMIT) {
    151             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    152         }
     149        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     150        if (SN > SN_LIMIT) {
     151            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     152        }
    153153    }
    154154    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
    155         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    156         if (SN > SN_LIMIT) {
    157             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    158         }
     155        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     156        if (SN > SN_LIMIT) {
     157            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     158        }
    159159    }
    160160
    … …  
    162162    // first go up:
    163163    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
    164         // mark the pixels in this row to the left, then the right
    165         for (int ix = 0; ix < pixels->numCols; ix++) {
    166             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    167             if (SN < SN_LIMIT) continue;
    168 
    169             bool valid = false;
    170             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
    171             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
    172             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
    173 
    174             if (!valid) continue;
    175             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    176         }
     164        // mark the pixels in this row to the left, then the right
     165        for (int ix = 0; ix < pixels->numCols; ix++) {
     166            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     167            if (SN < SN_LIMIT) continue;
     168
     169            bool valid = false;
     170            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
     171            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
     172            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
     173
     174            if (!valid) continue;
     175            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     176        }
    177177    }
    178178    // next go down:
    179179    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
    180         // mark the pixels in this row to the left, then the right
    181         for (int ix = 0; ix < pixels->numCols; ix++) {
    182             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    183             if (SN < SN_LIMIT) continue;
    184 
    185             bool valid = false;
    186             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
    187             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
    188             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
    189 
    190             if (!valid) continue;
    191             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    192         }
     180        // mark the pixels in this row to the left, then the right
     181        for (int ix = 0; ix < pixels->numCols; ix++) {
     182            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     183            if (SN < SN_LIMIT) continue;
     184
     185            bool valid = false;
     186            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
     187            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
     188            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
     189
     190            if (!valid) continue;
     191            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     192        }
    193193    }
    194194    return true;
    … …  
    201201    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
    202202    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    203    
     203
    204204    psImageMaskType maskVal = options->maskVal | options->markVal;
    205205
    … …  
    208208
    209209    for (int i = 0; i < sources->n; i++) {
    210         pmSource *source = sources->data[i];
    211         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
     210        pmSource *source = sources->data[i];
     211        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    212212
    213213        // replace object in image
    … …  
    216216        }
    217217
    218         // clear the mask bit and set the circular mask pixels
    219         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    220         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    221 
    222         // XXX can we test if psfMag is set and calculate only if needed?
    223         pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
    224        
    225         // clear the mask bit
    226         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     218        // clear the mask bit and set the circular mask pixels
     219        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     220        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
     221
     222        // XXX can we test if psfMag is set and calculate only if needed?
     223        pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
     224
     225        // clear the mask bit
     226        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    227227
    228228        // re-subtract the object, leave local sky
    229229        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    230230
    231         float apMag = -2.5*log10(source->moments->Sum);
    232         float dMag = source->psfMag - apMag;
    233        
    234         psVectorAppend (Ap, 100, dMag);
    235         psVectorAppend (ApErr, 100, source->errMag);
     231        float apMag = -2.5*log10(source->moments->Sum);
     232        float dMag = source->psfMag - apMag;
     233
     234        psVectorAppend (Ap, 100, dMag);
     235        psVectorAppend (ApErr, 100, source->errMag);
    236236    }
    237237
    … …  
    242242    psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
    243243    for (int i = 0; i < Ap->n; i++) {
    244         dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
     244        dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
    245245    }
    246246
    … …  
    266266    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");
    267267
    268     int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     268    int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS");
    269269    for (int i = 0; i < nRegions; i ++) {
    270         snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
    271         psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    272         psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
    273 
    274         psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
    275         psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
    276 
    277         // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
    278         psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
    279         psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
    280         psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
    281         psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
    282 
    283         if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
    284             psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    285             continue;
    286         }
    287        
    288         if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) {
    289             psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    290             continue;
    291         }
    292     }   
     270        snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
     271        psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, regionName);
     272        psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
     273
     274        psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
     275        psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
     276
     277        // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
     278        psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
     279        psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
     280        psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
     281        psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
     282
     283        if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
     284            psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     285            continue;
     286        }
     287
     288        if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) {
     289            psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     290            continue;
     291        }
     292    }
    293293
    294294    return true;
    … …  
    314314    for (psS32 i = 0 ; i < sources->n ; i++) {
    315315
    316         pmSource *source = (pmSource *) sources->data[i];
    317 
    318         // psfClumps are found for image subregions:
    319         // skip sources not in this region
    320         if (source->peak->x <  region->x0) continue;
    321         if (source->peak->x >= region->x1) continue;
    322         if (source->peak->y <  region->y0) continue;
    323         if (source->peak->y >= region->y1) continue;
    324 
    325         // skip source if it was already measured
    326         if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
    327             psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
    328             continue;
    329         }
    330 
    331         // source must have been subtracted
    332         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    333             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    334             psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
    335             continue;
    336         }
    337 
    338         // we are basically classifying by moments; use the default if not found
    339         psAssert (source->moments, "why is this source missing moments?");
    340         if (source->mode & noMoments) {
    341             Nskip ++;
    342             continue;
    343         }
    344 
    345         psF32 Mxx = source->moments->Mxx;
    346         psF32 Myy = source->moments->Myy;
     316        pmSource *source = (pmSource *) sources->data[i];
     317
     318        // psfClumps are found for image subregions:
     319        // skip sources not in this region
     320        if (source->peak->x <  region->x0) continue;
     321        if (source->peak->x >= region->x1) continue;
     322        if (source->peak->y <  region->y0) continue;
     323        if (source->peak->y >= region->y1) continue;
     324
     325        // skip source if it was already measured
     326        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     327            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     328            continue;
     329        }
     330
     331        // source must have been subtracted
     332        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     333            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
     334            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
     335            continue;
     336        }
     337
     338        // we are basically classifying by moments; use the default if not found
     339        psAssert (source->moments, "why is this source missing moments?");
     340        if (source->mode & noMoments) {
     341            Nskip ++;
     342            continue;
     343        }
     344
     345        psF32 Mxx = source->moments->Mxx;
     346        psF32 Myy = source->moments->Myy;
    347347
    348348        // replace object in image
    … …  
    351351        }
    352352
    353         // clear the mask bit and set the circular mask pixels
    354         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    355         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    356 
    357         // XXX can we test if psfMag is set and calculate only if needed?
    358         pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
    359 
    360         // clear the mask bit
    361         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     353        // clear the mask bit and set the circular mask pixels
     354        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     355        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
     356
     357        // XXX can we test if psfMag is set and calculate only if needed?
     358        pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
     359
     360        // clear the mask bit
     361        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    362362
    363363        // re-subtract the object, leave local sky
    364364        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    365365
    366         float apMag = -2.5*log10(source->moments->Sum);
    367         float dMag = source->psfMag - apMag;
    368         float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
    369 
    370         source->extNsigma = nSigma;
    371         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    372 
    373         // Anything within this region is a probably PSF-like object. Saturated stars may land
    374         // in this region, but are detected elsewhere on the basis of their peak value.
    375         bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
    376         if (isPSF) {
    377             Npsf ++;
    378             continue;
    379         }
    380 
    381         // Defects may not always match CRs from peak curvature analysis
    382         // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
    383         // XXX this rule is not great
    384         if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
    385             source->mode |= PM_SOURCE_MODE_DEFECT;
    386             Ncr ++;
    387             continue;
    388         }
    389 
    390         // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
    391         // just large saturated regions.
    392         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    393             Nsat ++;
    394             continue;
    395         }
    396 
    397         // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
    398         bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
    399         if (isEXT) {
    400             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    401             Next ++;
    402             continue;
    403         }
    404 
    405         psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
    406         Nmiss ++;
     366        float apMag = -2.5*log10(source->moments->Sum);
     367        float dMag = source->psfMag - apMag;
     368        float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
     369
     370        source->extNsigma = nSigma;
     371        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     372
     373        // Anything within this region is a probably PSF-like object. Saturated stars may land
     374        // in this region, but are detected elsewhere on the basis of their peak value.
     375        bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
     376        if (isPSF) {
     377            Npsf ++;
     378            continue;
     379        }
     380
     381        // Defects may not always match CRs from peak curvature analysis
     382        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
     383        // XXX this rule is not great
     384        if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
     385            source->mode |= PM_SOURCE_MODE_DEFECT;
     386            Ncr ++;
     387            continue;
     388        }
     389
     390        // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
     391        // just large saturated regions.
     392        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     393            Nsat ++;
     394            continue;
     395        }
     396
     397        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
     398        bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
     399        if (isEXT) {
     400            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     401            Next ++;
     402            continue;
     403        }
     404
     405        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
     406        Nmiss ++;
    407407    }
    408408
    … …  
    417417// no longer used by psphot.
    418418float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    419                         psImageMaskType maskVal, const pmModel *model, float Ro)
     419                        psImageMaskType maskVal, const pmModel *model, float Ro)
    420420{
    421421    psF32 *PAR = model->params->data.F32; // Model parameters
    … …  
    432432    float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
    433433    if (Q < 0.0) {
    434         // ellipse is imaginary
    435         return NAN;
     434        // ellipse is imaginary
     435        return NAN;
    436436    }
    437437
    … …  
    441441
    442442    for (int x = -radius; x <= radius; x++) {
    443         // Polynomial coefficients
    444         // XXX Should we be using the centre of the pixel as x or x+0.5?
    445         float A = PS_SQR (1.0 / syy);
    446         float B = x * sxy;
    447         float C = PS_SQR (x / sxx) - Ro;
    448         float T = PS_SQR(B) - 4*A*C;
    449         if (T < 0.0) {
    450             continue;
    451         }
    452 
    453         // y position in source frame
    454         float yP = (-B + sqrt (T)) / (2.0 * A);
    455         float yM = (-B - sqrt (T)) / (2.0 * A);
    456 
    457         // Get the closest pixel positions (image frame)
    458         int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
    459         int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    460         int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    461 
    462         if (xPix < 0 || xPix >= image->numCols) {
    463             continue;
    464         }
    465 
    466         if (yPixM >= 0 && yPixM < image->numRows &&
    467             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
    468             float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
    469             nSigma += dSigma;
    470             nPts++;
    471         }
    472 
    473         if (yPixM == yPixP) {
    474             continue;
    475         }
    476 
    477         if (yPixP >= 0 && yPixP < image->numRows &&
    478             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
    479             float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
    480             nSigma += dSigma;
    481             nPts++;
    482         }
     443        // Polynomial coefficients
     444        // XXX Should we be using the centre of the pixel as x or x+0.5?
     445        float A = PS_SQR (1.0 / syy);
     446        float B = x * sxy;
     447        float C = PS_SQR (x / sxx) - Ro;
     448        float T = PS_SQR(B) - 4*A*C;
     449        if (T < 0.0) {
     450            continue;
     451        }
     452
     453        // y position in source frame
     454        float yP = (-B + sqrt (T)) / (2.0 * A);
     455        float yM = (-B - sqrt (T)) / (2.0 * A);
     456
     457        // Get the closest pixel positions (image frame)
     458        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
     459        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     460        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     461
     462        if (xPix < 0 || xPix >= image->numCols) {
     463            continue;
     464        }
     465
     466        if (yPixM >= 0 && yPixM < image->numRows &&
     467            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
     468            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
     469            nSigma += dSigma;
     470            nPts++;
     471        }
     472
     473        if (yPixM == yPixP) {
     474            continue;
     475        }
     476
     477        if (yPixP >= 0 && yPixP < image->numRows &&
     478            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
     479            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
     480            nSigma += dSigma;
     481            nPts++;
     482        }
    483483    }
    484484    nSigma /= nPts;
    … …  
    491491    // XXX use an internal flag to mark sources which have already been measured
    492492    for (int i = 0; i < sources->n; i++) {
    493         pmSource *source = sources->data[i];
    494 
    495         // skip source if it was already measured
    496         if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
    497             psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
    498             continue;
    499         }
    500 
    501         // source must have been subtracted
    502         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    503             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    504             psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
    505             continue;
    506         }
    507 
    508         psF32 **resid  = source->pixels->data.F32;
    509         psF32 **variance = source->variance->data.F32;
    510         psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    511 
    512         // Integer position of peak
    513         int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
    514         int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
    515 
    516         // Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
    517         if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
    518             yPeak < 1 || yPeak > source->pixels->numRows - 2) {
    519             psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    520             continue;
    521         }
    522 
    523         // Skip sources with masked pixels.  These are mostly caught as DEFECT
    524         bool keep = true;
    525         for (int iy = -1; (iy <= +1) && keep; iy++) {
    526             for (int ix = -1; (ix <= +1) && keep; ix++) {
    527                 if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
    528                     keep = false;
    529                 }
    530             }
    531         }
    532         if (!keep) {
    533             psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
    534             continue;
    535         }
    536 
    537         // Compare the central pixel with those on either side, for the four possible lines through it.
    538 
    539         // Soften variances (add systematic error)
    540         float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances
    541 
    542         // Across the middle: y = 0
    543         float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
    544         float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
    545         float nX = cX / sqrtf(dcX + softening);
    546 
    547         // Up the centre: x = 0
    548         float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
    549         float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
    550         float nY = cY / sqrtf(dcY + softening);
    551 
    552         // Diagonal: x = y
    553         float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
    554         float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
    555         float nL = cL / sqrtf(dcL + softening);
    556 
    557         // Diagonal: x = - y
    558         float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
    559         float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
    560         float nR = cR / sqrtf(dcR + softening);
    561 
    562         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
    563         // Ndof = 4 ? (four measurements, no free parameters)
    564         // XXX this value is going to be biased low because of systematic errors.
    565         // we need to calibrate it somehow
    566         // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
    567 
    568         // not strictly accurate: overcounts the chisq contribution from the center pixel (by
    569         // factor of 4); also biases a bit low if any pixels are masked
    570         // XXX I am not sure I want to keep this value...
    571         source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
    572 
    573         float fCR = 0.0;
    574         int nCR = 0;
    575         if (nX > 0.0) {
    576             fCR += nX;
    577             nCR ++;
    578         }
    579         if (nY > 0.0) {
    580             fCR += nY;
    581             nCR ++;
    582         }
    583         if (nL > 0.0) {
    584             fCR += nL;
    585             nCR ++;
    586         }
    587         if (nR > 0.0) {
    588             fCR += nR;
    589             nCR ++;
    590         }
    591         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
    592         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    593 
    594         if (!isfinite(source->crNsigma)) {
    595             continue;
    596         }
    597 
    598         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    599         if (source->crNsigma > options->nSigmaCR) {
    600             source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    601             // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
    602             // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
    603         }
     493        pmSource *source = sources->data[i];
     494
     495        // skip source if it was already measured
     496        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     497            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     498            continue;
     499        }
     500
     501        // source must have been subtracted
     502        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     503            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
     504            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
     505            continue;
     506        }
     507
     508        psF32 **resid  = source->pixels->data.F32;
     509        psF32 **variance = source->variance->data.F32;
     510        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     511
     512        // Integer position of peak
     513        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
     514        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
     515
     516        // Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
     517        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
     518            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
     519            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
     520            continue;
     521        }
     522
     523        // Skip sources with masked pixels.  These are mostly caught as DEFECT
     524        bool keep = true;
     525        for (int iy = -1; (iy <= +1) && keep; iy++) {
     526            for (int ix = -1; (ix <= +1) && keep; ix++) {
     527                if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
     528                    keep = false;
     529                }
     530            }
     531        }
     532        if (!keep) {
     533            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
     534            continue;
     535        }
     536
     537        // Compare the central pixel with those on either side, for the four possible lines through it.
     538
     539        // Soften variances (add systematic error)
     540        float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances
     541
     542        // Across the middle: y = 0
     543        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
     544        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
     545        float nX = cX / sqrtf(dcX + softening);
     546
     547        // Up the centre: x = 0
     548        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
     549        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
     550        float nY = cY / sqrtf(dcY + softening);
     551
     552        // Diagonal: x = y
     553        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
     554        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
     555        float nL = cL / sqrtf(dcL + softening);
     556
     557        // Diagonal: x = - y
     558        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
     559        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
     560        float nR = cR / sqrtf(dcR + softening);
     561
     562        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
     563        // Ndof = 4 ? (four measurements, no free parameters)
     564        // XXX this value is going to be biased low because of systematic errors.
     565        // we need to calibrate it somehow
     566        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
     567
     568        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
     569        // factor of 4); also biases a bit low if any pixels are masked
     570        // XXX I am not sure I want to keep this value...
     571        source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
     572
     573        float fCR = 0.0;
     574        int nCR = 0;
     575        if (nX > 0.0) {
     576            fCR += nX;
     577            nCR ++;
     578        }
     579        if (nY > 0.0) {
     580            fCR += nY;
     581            nCR ++;
     582        }
     583        if (nL > 0.0) {
     584            fCR += nL;
     585            nCR ++;
     586        }
     587        if (nR > 0.0) {
     588            fCR += nR;
     589            nCR ++;
     590        }
     591        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
     592        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     593
     594        if (!isfinite(source->crNsigma)) {
     595            continue;
     596        }
     597
     598        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
     599        if (source->crNsigma > options->nSigmaCR) {
     600            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     601            // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
     602            // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
     603        }
    604604    }
    605605
    606606    // now that we have masked pixels associated with CRs, we can grow the mask
    607607    if (options->grow > 0) {
    608         bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
    609         psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
    610         psImageConvolveSetThreads(oldThreads);
    611         if (!newMask) {
    612             psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
    613             return false;
    614         }
    615         psFree(readout->mask);
    616         readout->mask = newMask;
     608        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
     609        psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
     610        psImageConvolveSetThreads(oldThreads);
     611        if (!newMask) {
     612            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
     613            return false;
     614        }
     615        psFree(readout->mask);
     616        readout->mask = newMask;
    617617    }
    618618    return true;

  • trunk/psphot/src/psphotSourceStats.c

    r25755 r25852  
    11# include "psphotInternal.h"
    22
    3 bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources);
     3bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    44
    55psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) {
    … …  
    6868
    6969    if (setWindow) {
    70         if (!psphotSetMomentsWindow(recipe, sources)) {
    71             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
    72             return NULL;
    73         }
     70        if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
     71            psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
     72            return NULL;
     73        }
    7474    }
    7575
    … …  
    242242
    243243// this function attempts to iteratively determine the best value for sigma of the moments weighting Gaussian
    244 bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources) {
     244bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources) {
    245245
    246246    bool status;
    … …  
    259259    for (int i = 0; i < 4; i++) {
    260260
    261         // XXX move max source number to config
    262         for (int j = 0; (j < sources->n) && (j < 400); j++) {
    263  
    264             pmSource *source = sources->data[j];
    265             psAssert (source->moments, "force moments to exist");
    266             source->moments->nPixels = 0;
    267 
    268             // skip faint sources for moments measurement
    269             if (source->peak->SN < MIN_SN) {
    270                 continue;
    271             }
    272 
    273             // measure basic source moments (no S/N clipping on input pixels)
    274             status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
    275         }
    276 
    277         // choose a grid scale that is a fixed fraction of the psf sigma^2
    278         psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
    279         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
    280         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
    281 
    282         // determine the PSF parameters from the source moment values
    283         pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
    284         psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
    285 
    286         psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
    287         psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    288         if (!regionMD) {
    289             regionMD = psMetadataAlloc();
    290             psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
    291             psFree (regionMD);
    292         }
    293         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
    294         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
    295         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
    296         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
    297            
    298         // psphotVisualPlotMoments (recipe, sources);
    299 
    300         Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     261        // XXX move max source number to config
     262        for (int j = 0; (j < sources->n) && (j < 400); j++) {
     263
     264            pmSource *source = sources->data[j];
     265            psAssert (source->moments, "force moments to exist");
     266            source->moments->nPixels = 0;
     267
     268            // skip faint sources for moments measurement
     269            if (source->peak->SN < MIN_SN) {
     270                continue;
     271            }
     272
     273            // measure basic source moments (no S/N clipping on input pixels)
     274            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
     275        }
     276
     277        // choose a grid scale that is a fixed fraction of the psf sigma^2
     278        psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
     279        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
     280        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
     281
     282        // determine the PSF parameters from the source moment values
     283        pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
     284        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
     285
     286        psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
     287        psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, "PSF.CLUMP.REGION.000");
     288        if (!regionMD) {
     289            regionMD = psMetadataAlloc();
     290            psMetadataAddMetadata (analysis, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
     291            psFree (regionMD);
     292        }
     293        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
     294        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
     295        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
     296        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
     297
     298        // psphotVisualPlotMoments (recipe, sources);
     299
     300        Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
    301301    }
    302302
    … …  
    307307    float maxS = Sout[0];
    308308    for (int i = 0; i < 4; i++) {
    309         minS = PS_MIN(Sout[i], minS);
    310         maxS = PS_MAX(Sout[i], maxS);
     309        minS = PS_MIN(Sout[i], minS);
     310        maxS = PS_MAX(Sout[i], maxS);
    311311    }
    312312    if (minS > 0.65) Sigma = sigma[3];
    … …  
    314314
    315315    for (int i = 0; i < 3; i++) {
    316         if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
    317         if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
    318         Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
     316        if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
     317        if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
     318        Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
    319319    }
    320320    psAssert (isfinite(Sigma), "did we miss a case?");
    321        
     321
    322322    // choose a grid scale that is a fixed fraction of the psf sigma^2
    323323    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));

  • trunk/psphot/src/psphotVisual.c

    r25755 r25852  
    2828    switch (channel) {
    2929      case 1:
    30         if (kapa1 == -1) {
    31             kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images");
    32             if (kapa1 == -1) {
    33                 fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    34                 pmVisualSetVisual(false);
    35             }
    36         }
    37         return kapa1;
     30        if (kapa1 == -1) {
     31            kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images");
     32            if (kapa1 == -1) {
     33                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     34                pmVisualSetVisual(false);
     35            }
     36        }
     37        return kapa1;
    3838      case 2:
    39         if (kapa2 == -1) {
    40             kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    41             if (kapa2 == -1) {
    42                 fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    43                 pmVisualSetVisual(false);
    44             }
    45         }
    46         return kapa2;
     39        if (kapa2 == -1) {
     40            kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots");
     41            if (kapa2 == -1) {
     42                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     43                pmVisualSetVisual(false);
     44            }
     45        }
     46        return kapa2;
    4747      case 3:
    48         if (kapa3 == -1) {
    49             kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps");
    50             if (kapa3 == -1) {
    51                 fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    52                 pmVisualSetVisual(false);
    53             }
    54         }
    55         return kapa3;
     48        if (kapa3 == -1) {
     49            kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps");
     50            if (kapa3 == -1) {
     51                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     52                pmVisualSetVisual(false);
     53            }
     54        }
     55        return kapa3;
    5656      default:
    57         psAbort ("unknown kapa channel");
     57        psAbort ("unknown kapa channel");
    5858    }
    5959    psAbort ("unknown kapa channel");
    … …  
    315315
    316316        // draw the top
    317         // XXX need to allow top (and bottom) to have more than one span
     317        // XXX need to allow top (and bottom) to have more than one span
    318318        span = footprint->spans->data[0];
    319319        overlay[Noverlay].type = KII_OVERLAY_LINE;
    … …  
    448448}
    449449
    450 bool psphotVisualPlotMoments (psMetadata *recipe, psArray *sources) {
     450bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources) {
    451451
    452452    bool status;
    … …  
    469469    float Ymin = 1000.0, Ymax = 0.0;
    470470    {
    471         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     471        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
    472472        for (int n = 0; n < nRegions; n++) {
    473473
    474474            char regionName[64];
    475475            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    476             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    477 
    478             float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     476            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     477
     478            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    479479            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    480480            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
    481481            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
    482482
    483             float X0 = psfX - 4.0*psfdX;
    484             float X1 = psfX + 4.0*psfdX;
    485             float Y0 = psfY - 4.0*psfdY;
    486             float Y1 = psfY + 4.0*psfdY;
    487 
    488             if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
    489             if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
    490             if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
    491             if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
    492         }
    493     }
    494     Xmin = PS_MAX(Xmin, -0.1); 
    495     Ymin = PS_MAX(Ymin, -0.1); 
     483            float X0 = psfX - 4.0*psfdX;
     484            float X1 = psfX + 4.0*psfdX;
     485            float Y0 = psfY - 4.0*psfdY;
     486            float Y1 = psfY + 4.0*psfdY;
     487
     488            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
     489            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     490            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
     491            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
     492        }
     493    }
     494    Xmin = PS_MAX(Xmin, -0.1);
     495    Ymin = PS_MAX(Ymin, -0.1);
    496496
    497497    // storage vectors for data to be plotted
    … …  
    653653        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
    654654
    655         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
    656         float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
     655        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
     656        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, analysis, "PSF_CLUMP_NSIGMA");
    657657
    658658        graphdata.color = KapaColorByName ("blue");
    659659        graphdata.style = 0;
    660660
    661         graphdata.xmin = Xmin;
    662         graphdata.ymin = Ymin;
    663         graphdata.xmax = Xmax;
    664         graphdata.ymax = Ymax;
    665         KapaSetLimits (myKapa, &graphdata);
     661        graphdata.xmin = Xmin;
     662        graphdata.ymin = Ymin;
     663        graphdata.xmax = Xmax;
     664        graphdata.ymax = Ymax;
     665        KapaSetLimits (myKapa, &graphdata);
    666666
    667667        for (int n = 0; n < nRegions; n++) {
    … …  
    669669            char regionName[64];
    670670            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    671             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     671            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
    672672
    673673            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    … …  
    926926            if (Xo == 0) {
    927927                // place source alone on this row
    928                 bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    929                 if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     928                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     929                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    930930                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    931931
    932                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     932                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    933933                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    934934
    935                 if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     935                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    936936
    937937                Yo += DY;
    … …  
    943943                Xo = 0;
    944944
    945                 bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    946                 if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     945                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     946                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    947947                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    948948
    949                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     949                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    950950                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    951951
    952                 if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     952                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    953953
    954954                Xo = DX;
    … …  
    957957        } else {
    958958            // extend this row
    959             bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    960             if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     959            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     960            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    961961            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    962962
    963             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     963            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    964964            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    965             if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     965            if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    966966
    967967            Xo += DX;
    … …  
    10201020        pmSource *source = sources->data[i];
    10211021
    1022         // only show "real" saturated stars (not defects)
     1022        // only show "real" saturated stars (not defects)
    10231023        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
    10241024        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
    … …  
    10651065        pmSource *source = sources->data[i];
    10661066
    1067         // only show "real" saturated stars (not defects)
     1067        // only show "real" saturated stars (not defects)
    10681068        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
    10691069        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
    1070         nSAT ++;
     1070        nSAT ++;
    10711071
    10721072        if (Xo + DX > NX) {
    … …  
    10741074            if (Xo == 0) {
    10751075                // place source alone on this row
    1076                 bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    1077                 if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     1076                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1077                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    10781078                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1079                 if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     1079                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    10801080
    10811081                Yo += DY;
    … …  
    10871087                Xo = 0;
    10881088
    1089                 bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    1090                 if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     1089                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1090                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    10911091                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1092                 if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     1092                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    10931093
    10941094                Xo = DX;
    … …  
    10971097        } else {
    10981098            // extend this row
    1099             bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    1100             if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     1099            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1100            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    11011101            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1102             if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     1102            if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    11031103
    11041104            Xo += DX;
    … …  
    14061406
    14071407        // if (source->type != type) continue;
    1408         if (mode) {
    1409             if (keep) {
    1410                 if (!(source->mode & mode)) continue;
    1411             } else {
    1412                 if (source->mode & mode) continue;
    1413             }
    1414         }
     1408        if (mode) {
     1409            if (keep) {
     1410                if (!(source->mode & mode)) continue;
     1411            } else {
     1412                if (source->mode & mode) continue;
     1413            }
     1414        }
    14151415
    14161416        pmMoments *moments = source->moments;
    … …  
    14631463}
    14641464
    1465 bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources) {
     1465bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources) {
    14661466
    14671467    bool status;
    … …  
    14821482    float Ymin = 1000.0, Ymax = 0.0;
    14831483    {
    1484         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     1484        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
    14851485        for (int n = 0; n < nRegions; n++) {
    14861486
    14871487            char regionName[64];
    14881488            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    1489             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    1490 
    1491             float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     1489            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     1490
     1491            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    14921492            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    14931493            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
    14941494            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
    14951495
    1496             float X0 = psfX - 10.0*psfdX;
    1497             float X1 = psfX + 10.0*psfdX;
    1498             float Y0 = psfY - 10.0*psfdY;
    1499             float Y1 = psfY + 10.0*psfdY;
    1500 
    1501             if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
    1502             if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
    1503             if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
    1504             if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
    1505         }
    1506     }
    1507     Xmin = PS_MAX(Xmin, -0.1); 
    1508     Ymin = PS_MAX(Ymin, -0.1); 
     1496            float X0 = psfX - 10.0*psfdX;
     1497            float X1 = psfX + 10.0*psfdX;
     1498            float Y0 = psfY - 10.0*psfdY;
     1499            float Y1 = psfY + 10.0*psfdY;
     1500
     1501            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
     1502            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     1503            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
     1504            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
     1505        }
     1506    }
     1507    Xmin = PS_MAX(Xmin, -0.1);
     1508    Ymin = PS_MAX(Ymin, -0.1);
    15091509
    15101510    // storage vectors for data to be plotted
    … …  
    15441544        if (source->moments == NULL) continue;
    15451545
    1546         if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
    1547             xCR->data.F32[nCR] = source->moments->Mxx;
    1548             yCR->data.F32[nCR] = source->moments->Myy;
    1549             mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum);
    1550             sCR->data.F32[nCR] = source->extNsigma;
    1551             nCR++;
    1552         }
    1553         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    1554             xSAT->data.F32[nSAT] = source->moments->Mxx;
    1555             ySAT->data.F32[nSAT] = source->moments->Myy;
    1556             mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum);
    1557             sSAT->data.F32[nSAT] = source->extNsigma;
    1558             nSAT++;
    1559         }
    1560         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    1561             xEXT->data.F32[nEXT] = source->moments->Mxx;
    1562             yEXT->data.F32[nEXT] = source->moments->Myy;
    1563             mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum);
    1564             sEXT->data.F32[nEXT] = source->extNsigma;
    1565             nEXT++;
    1566             continue;
    1567         }
    1568         if (source->mode & PM_SOURCE_MODE_DEFECT) {
    1569             xDEF->data.F32[nDEF] = source->moments->Mxx;
    1570             yDEF->data.F32[nDEF] = source->moments->Myy;
    1571             mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum);
    1572             sDEF->data.F32[nDEF] = source->extNsigma;
    1573             nDEF++;
    1574             continue;
    1575         }
    1576         if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
    1577             continue;
    1578         }
    1579         xPSF->data.F32[nPSF] = source->moments->Mxx;
    1580         yPSF->data.F32[nPSF] = source->moments->Myy;
    1581         mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum);
    1582         sPSF->data.F32[nPSF] = source->extNsigma;
    1583         nPSF++;
     1546        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
     1547            xCR->data.F32[nCR] = source->moments->Mxx;
     1548            yCR->data.F32[nCR] = source->moments->Myy;
     1549            mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum);
     1550            sCR->data.F32[nCR] = source->extNsigma;
     1551            nCR++;
     1552        }
     1553        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     1554            xSAT->data.F32[nSAT] = source->moments->Mxx;
     1555            ySAT->data.F32[nSAT] = source->moments->Myy;
     1556            mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum);
     1557            sSAT->data.F32[nSAT] = source->extNsigma;
     1558            nSAT++;
     1559        }
     1560        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     1561            xEXT->data.F32[nEXT] = source->moments->Mxx;
     1562            yEXT->data.F32[nEXT] = source->moments->Myy;
     1563            mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum);
     1564            sEXT->data.F32[nEXT] = source->extNsigma;
     1565            nEXT++;
     1566            continue;
     1567        }
     1568        if (source->mode & PM_SOURCE_MODE_DEFECT) {
     1569            xDEF->data.F32[nDEF] = source->moments->Mxx;
     1570            yDEF->data.F32[nDEF] = source->moments->Myy;
     1571            mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum);
     1572            sDEF->data.F32[nDEF] = source->extNsigma;
     1573            nDEF++;
     1574            continue;
     1575        }
     1576        if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
     1577            continue;
     1578        }
     1579        xPSF->data.F32[nPSF] = source->moments->Mxx;
     1580        yPSF->data.F32[nPSF] = source->moments->Myy;
     1581        mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum);
     1582        sPSF->data.F32[nPSF] = source->extNsigma;
     1583        nPSF++;
    15841584    }
    15851585    xSAT->n = nSAT;
    … …  
    18611861        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
    18621862
    1863         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
    1864         float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
     1863        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
     1864        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, analysis, "PSF_CLUMP_NSIGMA");
    18651865
    18661866        graphdata.color = KapaColorByName ("blue");
    18671867        graphdata.style = 0;
    18681868
    1869         graphdata.xmin = Xmin;
    1870         graphdata.ymin = Ymin;
    1871         graphdata.xmax = Xmax;
    1872         graphdata.ymax = Ymax;
    1873         KapaSetLimits (myKapa, &graphdata);
     1869        graphdata.xmin = Xmin;
     1870        graphdata.ymin = Ymin;
     1871        graphdata.xmax = Xmax;
     1872        graphdata.ymax = Ymax;
     1873        KapaSetLimits (myKapa, &graphdata);
    18741874
    18751875        for (int n = 0; n < nRegions; n++) {
    … …  
    18771877            char regionName[64];
    18781878            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    1879             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     1879            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
    18801880
    18811881            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    … …  
    20692069
    20702070    for (int i = 0; i < sources->n; i++) {
    2071         pmSource *source = sources->data[i];
    2072 
    2073         if (!source) continue;
    2074         if (!source->extpars) continue;
    2075         if (!source->extpars->profile) continue;
    2076         if (!source->extpars->petrosian_80) continue;
    2077 
    2078         pmSourceRadialProfile *profile = source->extpars->profile;
    2079         pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
    2080 
    2081         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    2082         overlay[Noverlay].x = source->peak->xf;
    2083         overlay[Noverlay].y = source->peak->yf;
    2084         overlay[Noverlay].dx = 2.0*petrosian->radius;
    2085         overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major;
    2086         overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
    2087         overlay[Noverlay].text = NULL;
    2088         Noverlay ++;
     2071        pmSource *source = sources->data[i];
     2072
     2073        if (!source) continue;
     2074        if (!source->extpars) continue;
     2075        if (!source->extpars->profile) continue;
     2076        if (!source->extpars->petrosian_80) continue;
     2077
     2078        pmSourceRadialProfile *profile = source->extpars->profile;
     2079        pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
     2080
     2081        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2082        overlay[Noverlay].x = source->peak->xf;
     2083        overlay[Noverlay].y = source->peak->yf;
     2084        overlay[Noverlay].dx = 2.0*petrosian->radius;
     2085        overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major;
     2086        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2087        overlay[Noverlay].text = NULL;
     2088        Noverlay ++;
    20892089        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    20902090
    2091         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    2092         overlay[Noverlay].x = source->peak->xf;
    2093         overlay[Noverlay].y = source->peak->yf;
    2094         overlay[Noverlay].dx = 4.0*petrosian->radius;
    2095         overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major;
    2096         overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
    2097         overlay[Noverlay].text = NULL;
    2098         Noverlay ++;
     2091        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2092        overlay[Noverlay].x = source->peak->xf;
     2093        overlay[Noverlay].y = source->peak->yf;
     2094        overlay[Noverlay].dx = 4.0*petrosian->radius;
     2095        overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major;
     2096        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2097        overlay[Noverlay].text = NULL;
     2098        Noverlay ++;
    20992099        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    21002100    }
Note: See TracChangeset for help on using the changeset viewer.