IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25359


Ignore:
Timestamp:
Sep 13, 2009, 4:56:34 PM (17 years ago)
Author:
eugene
Message:

use psfMag - momentsMag, momentsMxx, momentsMyy to classify source sizes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psphot/src/psphotSourceSize.c

    r25273 r25359  
    22# include <gsl/gsl_sf_gamma.h>
    33
    4 static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    5                                 psImageMaskType maskVal, const pmModel *model, float Ro);
     4bool psphotSourceSizePSF (float *ApResid, float *ApSysErr, psArray *sources, pmPSF *psf, psImageMaskType maskVal);
     5
     6bool psphotSourceSizeExtended (FILE *output, pmSource *source, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr);
     7
     8bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSyserr);
     9bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr);
     10
     11float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
     12                         psImageMaskType maskVal, const pmModel *model, float Ro);
    613
    714bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    … …  
    1421// deviation from the psf model at the r = FWHM/2 position
    1522
    16 bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first)
     23bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first)
    1724{
    1825    bool status;
    … …  
    3037    assert (status);
    3138
    32     float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
    33     assert (status);
     39    // float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
     40    // assert (status);
    3441
    3542    int grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
    … …  
    4552    }
    4653
    47     // loop over all source
     54    // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
     55    // of and object.  We need to model this distribution for the PSF stars before we can test
     56    // the significance for a specific object
     57    // XXX move this to the code that generates the PSF?
     58    // XXX store the results on pmPSF?
     59    float ApResid, ApSysErr;
     60    psphotSourceSizePSF (&ApResid, &ApSysErr, sources, psf, maskVal);
     61
     62    // classify the sources based on ApResid and Moments (extended sources)
     63    psphotSourceClass(readout, sources, recipe, psf, maskVal, ApResid, ApSysErr);
     64
     65    // classify the sources based on the CR test (place this in a function?)
    4866    for (int i = first; i < sources->n; i++) {
    49         pmSource *source = sources->data[i];
    50 
    51         // skip source if it was already measured
    52         if (isfinite(source->crNsigma)) {
    53             psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n");
    54             continue;
    55         }
    56 
    57         // source must have been subtracted
    58         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     67        pmSource *source = sources->data[i];
     68
     69        // skip source if it was already measured
     70        if (isfinite(source->crNsigma)) {
     71            psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n");
     72            continue;
     73        }
     74
     75        // source must have been subtracted
     76        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    5977            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    60             psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n");
    61             continue;
    62         }
    63 
    64         psF32 **resid  = source->pixels->data.F32;
    65         psF32 **variance = source->variance->data.F32;
    66         psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    67 
    68         // check for extendedness: measure the delta flux significance at the 1 sigma contour
    69         source->extNsigma = psphotModelContour(source->pixels, source->variance, source->maskObj, maskVal,
    70                                                source->modelPSF, 1.0);
    71 
    72         // XXX prevent a source from being both CR and EXT?
    73         if (source->extNsigma > EXT_NSIGMA_LIMIT) {
    74             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    75         }
    76 
    77         // Integer position of peak
    78         int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
    79         int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
    80 
    81         // XXX for now, skip sources which are too close to a boundary
    82         // XXX raise a flag?
    83         if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
    84             yPeak < 1 || yPeak > source->pixels->numRows - 2) {
     78            psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n");
     79            continue;
     80        }
     81
     82        psF32 **resid  = source->pixels->data.F32;
     83        psF32 **variance = source->variance->data.F32;
     84        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     85
     86        // Integer position of peak
     87        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
     88        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
     89
     90        // XXX for now, skip sources which are too close to a boundary
     91        // XXX raise a flag?
     92        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
     93            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
    8594            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    86             psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    87             continue;
    88         }
    89 
    90         // XXX for now, just skip any sources with masked pixels
    91         // XXX raise a flag?
    92         bool keep = true;
    93         for (int iy = -1; (iy <= +1) && keep; iy++) {
    94             for (int ix = -1; (ix <= +1) && keep; ix++) {
    95                 if (mask[yPeak+iy][xPeak+ix] & maskVal) {
    96                     keep = false;
    97                 }
    98             }
    99         }
    100         if (!keep) {
    101             psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
     95            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
     96            continue;
     97        }
     98
     99        // XXX for now, just skip any sources with masked pixels
     100        // XXX raise a flag?
     101        bool keep = true;
     102        for (int iy = -1; (iy <= +1) && keep; iy++) {
     103            for (int ix = -1; (ix <= +1) && keep; ix++) {
     104                if (mask[yPeak+iy][xPeak+ix] & maskVal) {
     105                    keep = false;
     106                }
     107            }
     108        }
     109        if (!keep) {
     110            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
    102111            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    103             continue;
    104         }
    105 
    106         // Compare the central pixel with those on either side, for the four possible lines through it.
    107 
    108         // Soften variances (add systematic error)
    109         float softening = soft * PS_SQR(source->peak->flux); // Softening for variances
    110 
    111         // Across the middle: y = 0
    112         float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
    113         float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
    114         float nX = cX / sqrtf(dcX + softening);
    115 
    116         // Up the centre: x = 0
    117         float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
    118         float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
    119         float nY = cY / sqrtf(dcY + softening);
    120 
    121         // Diagonal: x = y
    122         float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
    123         float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
    124         float nL = cL / sqrtf(dcL + softening);
    125 
    126         // Diagonal: x = - y
    127         float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
    128         float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
    129         float nR = cR / sqrtf(dcR + softening);
    130 
    131         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
    132         // Ndof = 4 ? (four measurements, no free parameters)
    133         // XXX this value is going to be biased low because of systematic errors.
    134         // we need to calibrate it somehow
    135         // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
    136 
    137         // not strictly accurate: overcounts the chisq contribution from the center pixel (by
    138         // factor of 4); also biases a bit low if any pixels are masked
    139         // XXX I am not sure I want to keep this value...
    140         source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
    141 
    142         float fCR = 0.0;
    143         int nCR = 0;
    144         if (nX > 0.0) {
    145             fCR += nX;
    146             nCR ++;
    147         }
    148         if (nY > 0.0) {
    149             fCR += nY;
    150             nCR ++;
    151         }
    152         if (nL > 0.0) {
    153             fCR += nL;
    154             nCR ++;
    155         }
    156         if (nR > 0.0) {
    157             fCR += nR;
    158             nCR ++;
    159         }
    160         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
    161         if (!isfinite(source->crNsigma)) {
    162             continue;
    163         }
    164 
    165         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    166         if (source->crNsigma > CR_NSIGMA_LIMIT) {
    167             // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
    168             // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
    169         }
    170     }
    171 
     112            continue;
     113        }
     114
     115        // Compare the central pixel with those on either side, for the four possible lines through it.
     116
     117        // Soften variances (add systematic error)
     118        float softening = soft * PS_SQR(source->peak->flux); // Softening for variances
     119
     120        // Across the middle: y = 0
     121        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
     122        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
     123        float nX = cX / sqrtf(dcX + softening);
     124
     125        // Up the centre: x = 0
     126        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
     127        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
     128        float nY = cY / sqrtf(dcY + softening);
     129
     130        // Diagonal: x = y
     131        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
     132        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
     133        float nL = cL / sqrtf(dcL + softening);
     134
     135        // Diagonal: x = - y
     136        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
     137        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
     138        float nR = cR / sqrtf(dcR + softening);
     139
     140        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
     141        // Ndof = 4 ? (four measurements, no free parameters)
     142        // XXX this value is going to be biased low because of systematic errors.
     143        // we need to calibrate it somehow
     144        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
     145
     146        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
     147        // factor of 4); also biases a bit low if any pixels are masked
     148        // XXX I am not sure I want to keep this value...
     149        source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
     150
     151        float fCR = 0.0;
     152        int nCR = 0;
     153        if (nX > 0.0) {
     154            fCR += nX;
     155            nCR ++;
     156        }
     157        if (nY > 0.0) {
     158            fCR += nY;
     159            nCR ++;
     160        }
     161        if (nL > 0.0) {
     162            fCR += nL;
     163            nCR ++;
     164        }
     165        if (nR > 0.0) {
     166            fCR += nR;
     167            nCR ++;
     168        }
     169        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
     170        if (!isfinite(source->crNsigma)) {
     171            continue;
     172        }
     173
     174        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
     175        if (source->crNsigma > CR_NSIGMA_LIMIT) {
     176            // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
     177            // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
     178        }
     179    }
     180
     181    // pause and wait for user input:
     182    // continue, save (provide name), ??
     183    char key[10];
     184    fprintf (stdout, "[c]ontinue? ");
     185    if (!fgets(key, 8, stdin)) {
     186        psWarning("Unable to read option");
     187    }
     188 
    172189    // now that we have masked pixels associated with CRs, we can grow the mask
    173190    if (grow > 0) {
    174         bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
    175         psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow);
    176         psImageConvolveSetThreads(oldThreads);
    177         if (!newMask) {
    178             psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
    179             return false;
    180         }
    181         psFree(readout->mask);
    182         readout->mask = newMask;
     191        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
     192        psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow);
     193        psImageConvolveSetThreads(oldThreads);
     194        if (!newMask) {
     195            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
     196            return false;
     197        }
     198        psFree(readout->mask);
     199        readout->mask = newMask;
    183200    }
    184201
    185202    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n",
    186               sources->n - first, psTimerMark ("psphot.size"));
     203              sources->n - first, psTimerMark ("psphot.size"));
    187204
    188205    psphotVisualPlotSourceSize (sources);
    … …  
    195212// deviation in sigmas.  This is measured on the residual image - should we ignore negative
    196213// deviations?
    197 static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    198                                 psImageMaskType maskVal, const pmModel *model, float Ro)
     214float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
     215                        psImageMaskType maskVal, const pmModel *model, float Ro)
    199216{
    200217    psF32 *PAR = model->params->data.F32; // Model parameters
    … …  
    211228    float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
    212229    if (Q < 0.0) {
    213         // ellipse is imaginary
    214         return NAN;
     230        // ellipse is imaginary
     231        return NAN;
    215232    }
    216233
    … …  
    220237
    221238    for (int x = -radius; x <= radius; x++) {
    222         // Polynomial coefficients
    223         // XXX Should we be using the centre of the pixel as x or x+0.5?
    224         float A = PS_SQR (1.0 / syy);
    225         float B = x * sxy;
    226         float C = PS_SQR (x / sxx) - Ro;
    227         float T = PS_SQR(B) - 4*A*C;
    228         if (T < 0.0) {
    229             continue;
    230         }
    231 
    232         // y position in source frame
    233         float yP = (-B + sqrt (T)) / (2.0 * A);
    234         float yM = (-B - sqrt (T)) / (2.0 * A);
    235 
    236         // Get the closest pixel positions (image frame)
    237         int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
    238         int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    239         int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    240 
    241         if (xPix < 0 || xPix >= image->numCols) {
    242             continue;
    243         }
    244 
    245         if (yPixM >= 0 && yPixM < image->numRows &&
    246             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
    247             float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
    248             nSigma += dSigma;
    249             nPts++;
    250         }
    251 
    252         if (yPixM == yPixP) {
    253             continue;
    254         }
    255 
    256         if (yPixP >= 0 && yPixP < image->numRows &&
    257             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
    258             float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
    259             nSigma += dSigma;
    260             nPts++;
    261         }
     239        // Polynomial coefficients
     240        // XXX Should we be using the centre of the pixel as x or x+0.5?
     241        float A = PS_SQR (1.0 / syy);
     242        float B = x * sxy;
     243        float C = PS_SQR (x / sxx) - Ro;
     244        float T = PS_SQR(B) - 4*A*C;
     245        if (T < 0.0) {
     246            continue;
     247        }
     248
     249        // y position in source frame
     250        float yP = (-B + sqrt (T)) / (2.0 * A);
     251        float yM = (-B - sqrt (T)) / (2.0 * A);
     252
     253        // Get the closest pixel positions (image frame)
     254        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
     255        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     256        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     257
     258        if (xPix < 0 || xPix >= image->numCols) {
     259            continue;
     260        }
     261
     262        if (yPixM >= 0 && yPixM < image->numRows &&
     263            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
     264            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
     265            nSigma += dSigma;
     266            nPts++;
     267        }
     268
     269        if (yPixM == yPixP) {
     270            continue;
     271        }
     272
     273        if (yPixP >= 0 && yPixP < image->numRows &&
     274            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
     275            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
     276            nSigma += dSigma;
     277            nPts++;
     278        }
    262279    }
    263280    nSigma /= nPts;
    … …  
    279296    pmFootprint *footprint = peak->footprint;
    280297    if (!footprint) {
    281         // if we have not footprint, use the old code to mask by isophot
    282         psphotMaskCosmicRay_Old (source, maskVal, crMask);
    283         return true;
     298        // if we have not footprint, use the old code to mask by isophot
     299        psphotMaskCosmicRay_Old (source, maskVal, crMask);
     300        return true;
    284301    }
    285302
    286303    if (!footprint->spans) {
    287         // if we have not footprint, use the old code to mask by isophot
    288         psphotMaskCosmicRay_Old (source, maskVal, crMask);
    289         return true;
     304        // if we have not footprint, use the old code to mask by isophot
     305        psphotMaskCosmicRay_Old (source, maskVal, crMask);
     306        return true;
    290307    }
    291308
    292309    // mask all of the pixels covered by the spans of the footprint
    293310    for (int j = 1; j < footprint->spans->n; j++) {
    294         pmSpan *span1 = footprint->spans->data[j];
    295 
    296         int iy = span1->y;
    297         int xs = span1->x0;
    298         int xe = span1->x1;
    299 
    300         for (int ix = xs; ix < xe; ix++) {
    301             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    302         }
     311        pmSpan *span1 = footprint->spans->data[j];
     312
     313        int iy = span1->y;
     314        int xs = span1->x0;
     315        int xe = span1->x1;
     316
     317        for (int ix = xs; ix < xe; ix++) {
     318            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     319        }
    303320    }
    304321    return true;
    … …  
    323340    // mark the pixels in this row to the left, then the right
    324341    for (int ix = xo; ix >= 0; ix--) {
    325         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    326         if (SN > SN_LIMIT) {
    327             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    328         }
     342        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     343        if (SN > SN_LIMIT) {
     344            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     345        }
    329346    }
    330347    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
    331         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    332         if (SN > SN_LIMIT) {
    333             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    334         }
     348        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     349        if (SN > SN_LIMIT) {
     350            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     351        }
    335352    }
    336353
    … …  
    338355    // first go up:
    339356    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
    340         // mark the pixels in this row to the left, then the right
    341         for (int ix = 0; ix < pixels->numCols; ix++) {
    342             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    343             if (SN < SN_LIMIT) continue;
    344 
    345             bool valid = false;
    346             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
    347             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
    348             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
    349 
    350             if (!valid) continue;
    351             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    352         }
     357        // mark the pixels in this row to the left, then the right
     358        for (int ix = 0; ix < pixels->numCols; ix++) {
     359            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     360            if (SN < SN_LIMIT) continue;
     361
     362            bool valid = false;
     363            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
     364            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
     365            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
     366
     367            if (!valid) continue;
     368            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     369        }
    353370    }
    354371    // next go down:
    355372    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
    356         // mark the pixels in this row to the left, then the right
    357         for (int ix = 0; ix < pixels->numCols; ix++) {
    358             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    359             if (SN < SN_LIMIT) continue;
    360 
    361             bool valid = false;
    362             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
    363             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
    364             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
    365 
    366             if (!valid) continue;
    367             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    368         }
    369     }
    370     return true;
    371 }
     373        // mark the pixels in this row to the left, then the right
     374        for (int ix = 0; ix < pixels->numCols; ix++) {
     375            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     376            if (SN < SN_LIMIT) continue;
     377
     378            bool valid = false;
     379            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
     380            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
     381            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
     382
     383            if (!valid) continue;
     384            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     385        }
     386    }
     387    return true;
     388}
     389
     390// is this source extended or not?
     391bool psphotSourceSizeExtended (FILE *output, pmSource *source, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) {
     392
     393    // XXX can we test if psfMag is set and calculate only if needed?
     394    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     395    pmSourceMagnitudes (source, psf, photMode, maskVal);
     396
     397    float apMag = -2.5*log10(source->moments->Sum);
     398    float dMag = source->psfMag - apMag;
     399   
     400    float nSigma = (dMag - ApResid) / hypot(source->errMag, ApSysErr);
     401
     402    fprintf (output, "%f %f : %f %f : %f %f : %f %f\n", source->peak->xf, source->peak->yf, nSigma, dMag, source->psfMag, source->errMag, source->moments->Mxx, source->moments->Myy);
     403
     404    return true;
     405}
     406
     407// model the apmifit distribution for the psf stars:
     408bool psphotSourceSizePSF (float *ApResid, float *ApSysErr, psArray *sources, pmPSF *psf, psImageMaskType maskVal) {
     409
     410    // select the psf stars
     411    psArray *psfstars = psArrayAllocEmpty (100);
     412
     413    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
     414    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
     415   
     416    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     417
     418    for (int i = 0; i < sources->n; i++) {
     419        pmSource *source = sources->data[i];
     420        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
     421        psArrayAdd (psfstars, 100, source);
     422
     423        // XXX can we test if psfMag is set and calculate only if needed?
     424        pmSourceMagnitudes (source, psf, photMode, maskVal);
     425       
     426        float apMag = -2.5*log10(source->moments->Sum);
     427        float dMag = source->psfMag - apMag;
     428       
     429        psVectorAppend (Ap, 100, dMag);
     430        psVectorAppend (ApErr, 100, source->errMag);
     431    }
     432
     433    // model the distribution as a mean or median value and a systematic error from that value:
     434    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     435    psVectorStats (stats, Ap, NULL, NULL, 0);
     436
     437    psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
     438    for (int i = 0; i < Ap->n; i++) {
     439        dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
     440    }
     441
     442    *ApResid = stats->robustMedian;
     443    *ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
     444    fprintf (stderr, "psf-sum: %f +/- %f\n", *ApResid, *ApSysErr);
     445
     446    return true;
     447}
     448
     449// classify sources based on the combination of psf-mag, Mxx, Myy
     450bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) {
     451
     452    bool status;
     453    pmPSFClump psfClump;
     454    char regionName[64];
     455
     456    int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     457    for (int i = 0; i < nRegions; i ++) {
     458        snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
     459        psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     460        psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
     461
     462        psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
     463        psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
     464
     465        // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
     466        psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
     467        psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
     468        psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
     469        psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
     470
     471        if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
     472            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);
     473            continue;
     474        }
     475       
     476        if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, maskVal, ApResid, ApSysErr)) {
     477            psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     478            continue;
     479        }
     480    }   
     481
     482    return true;
     483}
     484
     485bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) {
     486
     487    PS_ASSERT_PTR_NON_NULL(sources, false);
     488    PS_ASSERT_PTR_NON_NULL(recipe, false);
     489
     490    int Nsat  = 0;
     491    int Next  = 0;
     492    int Npsf  = 0;
     493    int Ncr   = 0;
     494    int Nmiss = 0;
     495    int Nskip = 0;
     496
     497    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
     498    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     499
     500    for (psS32 i = 0 ; i < sources->n ; i++) {
     501
     502        pmSource *source = (pmSource *) sources->data[i];
     503
     504        // psfClumps are found for image subregions:
     505        // skip sources not in this region
     506        if (source->peak->x <  region->x0) continue;
     507        if (source->peak->x >= region->x1) continue;
     508        if (source->peak->y <  region->y0) continue;
     509        if (source->peak->y >= region->y1) continue;
     510
     511        // we are basically classifying by moments; use the default if not found
     512        psAssert (source->moments, "why is this source missing moments?");
     513        if (source->mode & noMoments) {
     514            Nskip ++;
     515            continue;
     516        }
     517
     518        psF32 Mxx = source->moments->Mxx;
     519        psF32 Myy = source->moments->Myy;
     520
     521        // saturated star (determined in PSF fit)
     522        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     523            Nsat ++;
     524            continue;
     525        }
     526
     527        // XXX can we test if psfMag is set and calculate only if needed?
     528        pmSourceMagnitudes (source, psf, photMode, maskVal);
     529
     530        float apMag = -2.5*log10(source->moments->Sum);
     531        float dMag = source->psfMag - apMag;
     532        float nSigma = (dMag - ApResid) / hypot(source->errMag, ApSysErr);
     533
     534        source->extNsigma = nSigma;
     535
     536        // XXX these sigma cuts should be user-configured parameters
     537        bool isPSF = (fabs(nSigma) < 3.0) && (fabs(Mxx - psfClump->X) < 2.0*psfClump->dX) && (fabs(Myy - psfClump->Y) < 2.0*psfClump->dY);
     538        if (isPSF) {
     539            Npsf ++;
     540            continue;
     541        }
     542
     543        // XXX get the sign on nSigma right here and above
     544        bool isEXT = (nSigma < 3.0) && (Mxx > psfClump->X) && (Myy > psfClump->Y);
     545        if (isEXT) {
     546            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     547            Next ++;
     548            continue;
     549        }
     550
     551       
     552        // XXX possible or likely?  need to compare my standard CR test with this
     553        // Mark in both cases?
     554        if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
     555            Ncr ++;
     556            continue;
     557        }
     558
     559        Nmiss ++;
     560    }
     561
     562    psLogMsg("psModules.objects", PS_LOG_INFO, "Rough classifications: %d %d %d %d %d %d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
     563
     564    return true;
     565}
Note: See TracChangeset for help on using the changeset viewer.