IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 17, 2005, 10:41:29 AM (21 years ago)
Author:
magnier
Message:
  • moved the moments->SN calculation to pmSourceMoments
  • made more stringent requirement on min number of valid pixels for pmSourceMoments
  • dropped unneeded SN-related entries from pmSourceRoughClass
  • restricted test for SAT pixels to inner 5x5 square
  • added source->mode entries
  • modified source->type entries to match new options
  • added the PM_SOURCE_FITTED flag
  • dropped the weight modification based on model deviation in pmSourceFitModel (was introducing noise, not rejecting it)
File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_rel9_b0/psModules/src/objects/pmObjects.c

    r5727 r5806  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-07 20:01:05 $
     8 *  @version $Revision: 1.3.4.2 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-12-17 20:41:29 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    134134    tmp->modelPSF = NULL;
    135135    tmp->modelFLT = NULL;
    136     tmp->type = 0;
     136    tmp->type = PM_SOURCE_UNKNOWN;
     137    tmp->mode = PM_SOURCE_DEFAULT;
    137138    psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
    138139
    … …  
    749750    psS32 numPixels = 0;
    750751    psF32 Sum = 0.0;
     752    psF32 Var = 0.0;
    751753    psF32 X1 = 0.0;
    752754    psF32 Y1 = 0.0;
    … …  
    786788
    787789            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     790            psF32 wDiff = source->weight->data.F32[row][col];
    788791
    789792            // XXX EAM : check for valid S/N in pixel
    … …  
    793796            }
    794797
     798            Var += wDiff;
    795799            Sum += pDiff;
    796800            X1  += xDiff * pDiff;
    … …  
    805809        }
    806810    }
     811
     812    // if we have less than (1/4) of the possible pixels, force a retry
    807813    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    808     if ((numPixels < 3) || (Sum <= 0)) {
    809         psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
     814    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
     815        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
    810816        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    811817        return (false);
    … …  
    824830    y = Y1/Sum;
    825831    if ((fabs(x) > radius) || (fabs(y) > radius)) {
    826         psTrace (".psModules.pmSourceMoments", 5,
     832        psTrace (".psModules.pmSourceMoments", 3,
    827833                 "large centroid swing; invalid peak %d, %d\n",
    828834                 source->peak->x, source->peak->y);
    … …  
    837843    source->moments->Sxy = XY/Sum - x*y;
    838844    source->moments->Sum = Sum;
     845    source->moments->SN  = Sum / sqrt(Var);
    839846    source->moments->Peak = peakPixel;
    840847    source->moments->nPixels = numPixels;
    … …  
    10641071XXX: How can this function ever return FALSE?
    10651072 
    1066 XXX EAM : add the saturated mask value to metadata
     1073EAM: I moved S/N calculation to pmSourceMoments, using weight image
    10671074*****************************************************************************/
    10681075
    … …  
    10791086    int Ncr      = 0;
    10801087    int Nsatstar = 0;
    1081     psRegion allArray = psRegionSet (0, 0, 0, 0);
    1082 
     1088    // psRegion allArray = psRegionSet (0, 0, 0, 0);
     1089    psRegion inner;
     1090
     1091    // report stats on S/N values for star-like objects
    10831092    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    10841093    starsn->n = 0;
    … …  
    10861095    // check return status value (do these exist?)
    10871096    bool status;
    1088     psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
    1089     psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
    10901097    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
    1091     // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
    10921098
    10931099    // XXX allow clump size to be scaled relative to sigmas?
    … …  
    11021108        psF32 sigY = tmpSrc->moments->Sy;
    11031109
    1104         // calculate and save signal-to-noise estimates
    1105         psF32 S  = tmpSrc->moments->Sum;
    1106         psF32 A  = 4 * M_PI * sigX * sigY;
    1107         psF32 B  = tmpSrc->moments->Sky;
    1108         psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
    1109         psF32 SN = (S * sqrt(GAIN) / RT);
    1110         tmpSrc->moments->SN = SN;
    1111 
    11121110        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1113         int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
     1111        inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
     1112        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
    11141113
    11151114        // saturated star (size consistent with PSF or larger)
    … …  
    11171116        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
    11181117        if ((Nsatpix > 1) && big) {
    1119             tmpSrc->type = PM_SOURCE_SATSTAR;
     1118            tmpSrc->type = PM_SOURCE_STAR;
     1119            tmpSrc->mode = PM_SOURCE_SATSTAR;
    11201120            Nsatstar ++;
    11211121            continue;
    … …  
    11251125        if (Nsatpix > 1) {
    11261126            tmpSrc->type = PM_SOURCE_SATURATED;
     1127            tmpSrc->mode = PM_SOURCE_DEFAULT;
    11271128            Nsat ++;
    11281129            continue;
    … …  
    11341135        if ((sigX < 0.05) || (sigY < 0.05)) {
    11351136            tmpSrc->type = PM_SOURCE_DEFECT;
     1137            tmpSrc->mode = PM_SOURCE_DEFAULT;
    11361138            Ncr ++;
    11371139            continue;
    … …  
    11411143        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    11421144            tmpSrc->type = PM_SOURCE_GALAXY;
     1145            tmpSrc->mode = PM_SOURCE_DEFAULT;
    11431146            Ngal ++;
    11441147            continue;
    … …  
    11461149
    11471150        // the rest are probable stellar objects
    1148         starsn->data.F32[starsn->n] = SN;
     1151        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
    11491152        starsn->n ++;
    11501153        Nstar ++;
    11511154
    1152         // PSF star (within 1.5 sigma of clump center
     1155        // PSF star (within 1.5 sigma of clump center, S/N > limit)
    11531156        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    1154         if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
    1155             tmpSrc->type = PM_SOURCE_PSFSTAR;
     1157        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
     1158            tmpSrc->type = PM_SOURCE_STAR;
     1159            tmpSrc->mode = PM_SOURCE_PSFSTAR;
    11561160            Npsf ++;
    11571161            continue;
    … …  
    11591163
    11601164        // random type of star
    1161         tmpSrc->type = PM_SOURCE_OTHER;
     1165        tmpSrc->type = PM_SOURCE_STAR;
     1166        tmpSrc->mode = PM_SOURCE_DEFAULT;
    11621167    }
    11631168
    … …  
    16981703    }
    16991704
     1705    source->mode |= PM_SOURCE_FITTED;
     1706
    17001707    psFree(x);
    17011708    psFree(y);
    … …  
    17291736    psBool onPic     = true;
    17301737    psBool rc        = true;
    1731     psF32  Ro, ymodel;
     1738    // psF32  Ro, ymodel;
    17321739
    17331740    psVector *params = model->params;
    … …  
    17391746    // XXX EAM : I need to use the sky value to constrain the weight model
    17401747    int nParams = PSF ? params->n - 4 : params->n;
    1741     psF32 So = params->data.F32[0];
     1748    // psF32 So = params->data.F32[0];
    17421749
    17431750    // find the number of valid pixels
    … …  
    17761783
    17771784                // compare observed flux to model flux to adjust weight
    1778                 ymodel = modelFunc (NULL, model->params, coord);
     1785                // ymodel = modelFunc (NULL, model->params, coord);
    17791786
    17801787                // this test enhances the weight based on deviation from the model flux
    1781                 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
     1788                // Ro = 1.0 + sqrt (fabs (y->data.F32[tmpCnt] - ymodel) / (PS_SQR(ymodel - So) + 1.0));
    17821789
    17831790                // psMinimizeLMChi2_EAM takes wt = 1/dY^2
    … …  
    17851792                    yErr->data.F32[tmpCnt] = 0.0;
    17861793                } else {
    1787                     yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
     1794                    // yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
     1795                    yErr->data.F32[tmpCnt] = 1.0 / source->weight->data.F32[i][j];
    17881796                }
    17891797                tmpCnt++;
    … …  
    18651873        model->status = PM_MODEL_OFFIMAGE;
    18661874    }
     1875
     1876    source->mode |= PM_SOURCE_FITTED;
    18671877
    18681878    psFree(x);
Note: See TracChangeset for help on using the changeset viewer.