IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 5, 2010, 1:38:43 PM (16 years ago)
Author:
eugene
Message:

updates to psphot APIs to enable stack photometry

Location:
branches/eam_branches/20091201/psphot
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psphot

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

    r26646 r26788  
    2020bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
    2121bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
    22 bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    23 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    24 bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
    25 
    26 bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
    27 
    28     for (int i = 0; i < sources->n; i++) {
    29         pmSource *source = sources->data[i];
    30         pmPeak *peak = source->peak;
    31         pmFootprint *footprint = peak->footprint;
    32         if (!footprint) continue;
    33         for (int j = 0; j < footprint->spans->n; j++) {
    34             pmSpan *sp = footprint->spans->data[j];
    35             psAssert (sp, "missing span");
    36         }
    37     }
    38     return true;
    39 }
     22bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
     23bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
    4024
    4125// we need to call this function after sources have been fitted to the PSF model and
    … …  
    4529// deviation from the psf model at the r = FWHM/2 position
    4630
    47 // XXX use an internal flag to mark sources which have already been measured
    48 bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first)
     31// for now, let's store the detections on the readout->analysis for each readout
     32bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize)
     33{
     34    bool status = true;
     35
     36    // select the appropriate recipe information
     37    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     38    psAssert (recipe, "missing recipe?");
     39
     40    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     41    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     42
     43    // loop over the available readouts
     44    for (int i = 0; i < num; i++) {
     45        if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
     46            psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i);
     47            return false;
     48        }
     49    }
     50    return true;
     51}
     52
     53// this function use an internal flag to mark sources which have already been measured
     54bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize)
    4955{
    5056    bool status;
    … …  
    5258
    5359    psTimerStart ("psphot.size");
     60
     61    // find the currently selected readout
     62    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     63    psAssert (file, "missing file?");
     64
     65    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     66    psAssert (readout, "missing readout?");
     67
     68    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     69    psAssert (detections, "missing detections?");
     70
     71    psArray *sources = detections->allSources;
     72    psAssert (sources, "missing sources?");
     73
     74    if (!sources->n) {
     75        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     76        return true;
     77    }
     78
     79    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     80    psAssert (psf, "missing psf?");
    5481
    5582    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    91118    // XXX move this to the code that generates the PSF?
    92119    // XXX store the results on pmPSF?
    93     psphotSourceSizePSF (&options, sources, psf);
     120
     121    // XXX this should only be done on the first pass (ie, if we have newSources or allSources?)
     122    if (getPSFsize) {
     123        psphotSourceSizePSF (&options, sources, psf);
     124    }
    94125
    95126    // classify the sources based on ApResid and Moments (extended sources)
     127    // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    96128    psphotSourceClass(readout, sources, recipe, psf, &options);
    97129
     130    // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    98131    psphotSourceSizeCR (readout, sources, &options);
    99132
    100     psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size"));
     133    // XXX fix this (was source->n  - first)
     134    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
    101135
    102136    psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    103137    psphotVisualShowSourceSize (readout, sources);
    104138    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
    105 
    106     return true;
    107 }
    108 
    109 // This attempt to mask the cosmic rays used the isophotal boundary
    110 bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    111 
    112     // replace the source flux
    113     pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    114 
    115     // flag this as a CR
    116     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    117     pmPeak *peak = source->peak;
    118     psAssert (peak, "NULL peak");
    119 
    120     // grab the matching footprint
    121     pmFootprint *footprint = peak->footprint;
    122     if (!footprint) {
    123       psTrace("psphot.czw",2,"Using isophot CR mask code.");
    124      
    125         // if we have not footprint, use the old code to mask by isophot
    126         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    127         return true;
    128     }
    129 
    130     if (!footprint->spans) {
    131       psTrace("psphot.czw",2,"Using isophot CR mask code.");
    132      
    133         // if we have no footprint, use the old code to mask by isophot
    134         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    135         return true;
    136     }
    137     psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    138     // mask all of the pixels covered by the spans of the footprint
    139     for (int j = 1; j < footprint->spans->n; j++) {
    140         pmSpan *span1 = footprint->spans->data[j];
    141 
    142         int iy = span1->y;
    143         int xs = span1->x0;
    144         int xe = span1->x1;
    145 
    146         for (int ix = xs; ix < xe; ix++) {
    147             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    148         }
    149     }
    150     return true;
    151 }
    152 
    153 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    154 
    155     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    156     pmPeak *peak = source->peak;
    157     psAssert (peak, "NULL peak");
    158 
    159     psImage *mask   = source->maskView;
    160     psImage *pixels = source->pixels;
    161     psImage *variance = source->variance;
    162 
    163     // XXX This should be a recipe variable
    164 # define SN_LIMIT 5.0
    165 
    166     int xo = peak->x - pixels->col0;
    167     int yo = peak->y - pixels->row0;
    168 
    169     // mark the pixels in this row to the left, then the right
    170     for (int ix = xo; ix >= 0; ix--) {
    171         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    172         if (SN > SN_LIMIT) {
    173             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    174         }
    175     }
    176     for (int ix = xo + 1; ix < pixels->numCols; ix++) {
    177         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    178         if (SN > SN_LIMIT) {
    179             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    180         }
    181     }
    182 
    183     // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
    184     // first go up:
    185     for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
    186         // mark the pixels in this row to the left, then the right
    187         for (int ix = 0; ix < pixels->numCols; ix++) {
    188             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    189             if (SN < SN_LIMIT) continue;
    190 
    191             bool valid = false;
    192             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
    193             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
    194             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
    195 
    196             if (!valid) continue;
    197             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    198         }
    199     }
    200     // next go down:
    201     for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
    202         // mark the pixels in this row to the left, then the right
    203         for (int ix = 0; ix < pixels->numCols; ix++) {
    204             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    205             if (SN < SN_LIMIT) continue;
    206 
    207             bool valid = false;
    208             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
    209             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
    210             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
    211 
    212             if (!valid) continue;
    213             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    214         }
    215     }
     139    psphotVisualShowSatStars (recipe, psf, sources);
     140
    216141    return true;
    217142}
    … …  
    254179        float dMag = source->psfMag - apMag;
    255180
    256         psVectorAppend (Ap, 100, dMag);
    257         psVectorAppend (ApErr, 100, source->errMag);
     181        psVectorAppend (Ap, dMag);
     182        psVectorAppend (ApErr, source->errMag);
    258183    }
    259184
    … …  
    454379}
    455380
    456 // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
    457 // deviation in sigmas.  This is measured on the residual image - should we ignore negative
    458 // deviations?  NOTE: This function was an early attempt to classify extended objects, and is
    459 // no longer used by psphot.
    460 float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    461                          psImageMaskType maskVal, const pmModel *model, float Ro)
    462 {
    463     psF32 *PAR = model->params->data.F32; // Model parameters
    464     float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters
    465 
    466     // We treat the contour as an ellipse:
    467     // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY
    468     // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
    469     // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro
    470     // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so:
    471     // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]
    472 
    473     // min/max value of x is where B^2 - 4AC = 0; solve this for x
    474     float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
    475     if (Q < 0.0) {
    476         // ellipse is imaginary
    477         return NAN;
    478     }
    479 
    480     int radius = sqrtf(Q) + 0.5;        // Radius of ellipse
    481     int nPts = 0;                       // Number of points in ellipse
    482     float nSigma = 0.0;                 //
    483 
    484     for (int x = -radius; x <= radius; x++) {
    485         // Polynomial coefficients
    486         // XXX Should we be using the centre of the pixel as x or x+0.5?
    487         float A = PS_SQR (1.0 / syy);
    488         float B = x * sxy;
    489         float C = PS_SQR (x / sxx) - Ro;
    490         float T = PS_SQR(B) - 4*A*C;
    491         if (T < 0.0) {
    492             continue;
    493         }
    494 
    495         // y position in source frame
    496         float yP = (-B + sqrt (T)) / (2.0 * A);
    497         float yM = (-B - sqrt (T)) / (2.0 * A);
    498 
    499         // Get the closest pixel positions (image frame)
    500         int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
    501         int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    502         int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    503 
    504         if (xPix < 0 || xPix >= image->numCols) {
    505             continue;
    506         }
    507 
    508         if (yPixM >= 0 && yPixM < image->numRows &&
    509             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
    510             float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
    511             nSigma += dSigma;
    512             nPts++;
    513         }
    514 
    515         if (yPixM == yPixP) {
    516             continue;
    517         }
    518 
    519         if (yPixP >= 0 && yPixP < image->numRows &&
    520             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
    521             float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
    522             nSigma += dSigma;
    523             nPts++;
    524         }
    525     }
    526     nSigma /= nPts;
    527     return nSigma;
    528 }
    529 
    530381// given an object suspected to be a defect, generate a pixel mask using the Lapacian transform
    531382// if enough of the object is detected as 'sharp', consider the object a cosmic ray
    … …  
    534385    for (int i = 0; i < sources->n; i++) {
    535386        pmSource *source = sources->data[i];
     387
     388        // skip source if it was already measured
     389        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     390            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     391            continue;
     392        }
    536393
    537394        // Integer position of peak
    … …  
    582439// does no repair or recovery of the CR pixels, it only masks them out.  My test code can be
    583440// found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c
    584 bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
     441bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
    585442
    586443    // Get the actual images and information about the peak.
    … …  
    760617}
    761618
     619bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
     620
     621    for (int i = 0; i < sources->n; i++) {
     622        pmSource *source = sources->data[i];
     623        pmPeak *peak = source->peak;
     624        pmFootprint *footprint = peak->footprint;
     625        if (!footprint) continue;
     626        for (int j = 0; j < footprint->spans->n; j++) {
     627            pmSpan *sp = footprint->spans->data[j];
     628            psAssert (sp, "missing span");
     629        }
     630    }
     631    return true;
     632}
     633
     634/**** ------ old versions of cosmic ray masking ----- ****/
     635
     636bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
     637
     638// This attempt to mask the cosmic rays used the isophotal boundary
     639bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
     640
     641    // replace the source flux
     642    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     643
     644    // flag this as a CR
     645    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     646    pmPeak *peak = source->peak;
     647    psAssert (peak, "NULL peak");
     648
     649    // grab the matching footprint
     650    pmFootprint *footprint = peak->footprint;
     651    if (!footprint) {
     652      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     653     
     654        // if we have not footprint, use the old code to mask by isophot
     655        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     656        return true;
     657    }
     658
     659    if (!footprint->spans) {
     660      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     661     
     662        // if we have no footprint, use the old code to mask by isophot
     663        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     664        return true;
     665    }
     666    psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     667    // mask all of the pixels covered by the spans of the footprint
     668    for (int j = 1; j < footprint->spans->n; j++) {
     669        pmSpan *span1 = footprint->spans->data[j];
     670
     671        int iy = span1->y;
     672        int xs = span1->x0;
     673        int xe = span1->x1;
     674
     675        for (int ix = xs; ix < xe; ix++) {
     676            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     677        }
     678    }
     679    return true;
     680}
     681
     682bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
     683
     684    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     685    pmPeak *peak = source->peak;
     686    psAssert (peak, "NULL peak");
     687
     688    psImage *mask   = source->maskView;
     689    psImage *pixels = source->pixels;
     690    psImage *variance = source->variance;
     691
     692    // XXX This should be a recipe variable
     693# define SN_LIMIT 5.0
     694
     695    int xo = peak->x - pixels->col0;
     696    int yo = peak->y - pixels->row0;
     697
     698    // mark the pixels in this row to the left, then the right
     699    for (int ix = xo; ix >= 0; ix--) {
     700        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     701        if (SN > SN_LIMIT) {
     702            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     703        }
     704    }
     705    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
     706        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     707        if (SN > SN_LIMIT) {
     708            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     709        }
     710    }
     711
     712    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
     713    // first go up:
     714    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
     715        // mark the pixels in this row to the left, then the right
     716        for (int ix = 0; ix < pixels->numCols; ix++) {
     717            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     718            if (SN < SN_LIMIT) continue;
     719
     720            bool valid = false;
     721            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
     722            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
     723            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
     724
     725            if (!valid) continue;
     726            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     727        }
     728    }
     729    // next go down:
     730    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
     731        // mark the pixels in this row to the left, then the right
     732        for (int ix = 0; ix < pixels->numCols; ix++) {
     733            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     734            if (SN < SN_LIMIT) continue;
     735
     736            bool valid = false;
     737            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
     738            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
     739            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
     740
     741            if (!valid) continue;
     742            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     743        }
     744    }
     745    return true;
     746}
     747
     748// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
     749// deviation in sigmas.  This is measured on the residual image - should we ignore negative
     750// deviations?  NOTE: This function was an early attempt to classify extended objects, and is
     751// no longer used by psphot.
     752float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
     753                         psImageMaskType maskVal, const pmModel *model, float Ro)
     754{
     755    psF32 *PAR = model->params->data.F32; // Model parameters
     756    float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters
     757
     758    // We treat the contour as an ellipse:
     759    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY
     760    // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
     761    // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro
     762    // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so:
     763    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]
     764
     765    // min/max value of x is where B^2 - 4AC = 0; solve this for x
     766    float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
     767    if (Q < 0.0) {
     768        // ellipse is imaginary
     769        return NAN;
     770    }
     771
     772    int radius = sqrtf(Q) + 0.5;        // Radius of ellipse
     773    int nPts = 0;                       // Number of points in ellipse
     774    float nSigma = 0.0;                 //
     775
     776    for (int x = -radius; x <= radius; x++) {
     777        // Polynomial coefficients
     778        // XXX Should we be using the centre of the pixel as x or x+0.5?
     779        float A = PS_SQR (1.0 / syy);
     780        float B = x * sxy;
     781        float C = PS_SQR (x / sxx) - Ro;
     782        float T = PS_SQR(B) - 4*A*C;
     783        if (T < 0.0) {
     784            continue;
     785        }
     786
     787        // y position in source frame
     788        float yP = (-B + sqrt (T)) / (2.0 * A);
     789        float yM = (-B - sqrt (T)) / (2.0 * A);
     790
     791        // Get the closest pixel positions (image frame)
     792        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
     793        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     794        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     795
     796        if (xPix < 0 || xPix >= image->numCols) {
     797            continue;
     798        }
     799
     800        if (yPixM >= 0 && yPixM < image->numRows &&
     801            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
     802            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
     803            nSigma += dSigma;
     804            nPts++;
     805        }
     806
     807        if (yPixM == yPixP) {
     808            continue;
     809        }
     810
     811        if (yPixP >= 0 && yPixP < image->numRows &&
     812            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
     813            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
     814            nSigma += dSigma;
     815            nPts++;
     816        }
     817    }
     818    nSigma /= nPts;
     819    return nSigma;
     820}
     821
    762822// this was an old attempt to identify cosmic rays based on the peak curvature
    763823bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
    … …  
    893953    return true;
    894954}
     955
Note: See TracChangeset for help on using the changeset viewer.