IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10049


Ignore:
Timestamp:
Nov 17, 2006, 1:01:30 PM (19 years ago)
Author:
magnier
Message:

cleanups; added moment recalculation for satstars with larger box

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmSource.c

    r9880 r10049  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-11-07 09:06:21 $
     8 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-11-17 23:01:30 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3131static void sourceFree(pmSource *tmp)
    3232{
    33     psTrace("psModules.objects", 4, "---- %s() begin ----\n", __func__);
     33    psTrace(__func__, 5, "---- begin ----\n");
    3434    psFree(tmp->peak);
    3535    psFree(tmp->pixels);
    … …  
    4040    psFree(tmp->modelEXT);
    4141    psFree(tmp->blends);
    42     psTrace("psModules.objects", 4, "---- %s() end ----\n", __func__);
     42    psTrace(__func__, 5, "---- end ----\n");
    4343}
    4444
    … …  
    6565pmSource *pmSourceAlloc()
    6666{
    67     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
     67    psTrace(__func__, 5, "---- begin ----\n");
    6868    static int id = 1;
    6969    pmSource *source = (pmSource *) psAlloc(sizeof(pmSource));
    … …  
    9090    source->pixWeight = NAN;
    9191
    92     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
     92    psTrace(__func__, 5, "---- end ----\n");
    9393    return(source);
    9494}
    … …  
    198198
    199199/******************************************************************************
    200     pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
     200    pmSourcePSFClump(source, recipe): Find the likely PSF clump in the
    201201    sigma-x, sigma-y plane. return 0,0 clump in case of error.
    202202*****************************************************************************/
    203203
    204204// XXX EAM include a S/N cutoff in selecting the sources?
    205 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
    206 {
    207     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
     205pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *recipe)
     206{
     207    psTrace(__func__, 5, "---- begin ----\n");
    208208
    209209    # define NPIX 10
    … …  
    215215
    216216    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
    217     PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
     217    PS_ASSERT_PTR_NON_NULL(recipe, emptyClump);
    218218
    219219    // find the sigmaX, sigmaY clump
    … …  
    224224        bool status;
    225225
    226         psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
     226        psF32 SX_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_SX_MAX");
    227227        if (!status)
    228228            SX_MAX = 10.0;
    229         psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
     229        psF32 SY_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_SY_MAX");
    230230        if (!status)
    231231            SY_MAX = 10.0;
    232         psF32 AR_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_AR_MAX");
     232        psF32 AR_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_AR_MAX");
    233233        if (!status)
    234234            AR_MAX =  3.0;
    … …  
    291291        stats = psImageStats (stats, splane, NULL, 0);
    292292        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
    293         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
    294 
    295         const bool keep_psf_clump = psMetadataLookupBool(NULL, metadata, "KEEP_PSF_CLUMP");
     293        psTrace (__func__, 2, "clump threshold is %f\n", stats[0].max/2);
     294
     295        const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
    296296        if (keep_psf_clump)
    297297        {
    298             psMetadataAdd(metadata, PS_LIST_TAIL,
     298            psMetadataAdd(recipe, PS_LIST_TAIL,
    299299                          "PSF_CLUMP", PS_DATA_IMAGE, "Image of PSF coefficients", splane);
    300300        }
    … …  
    323323        psArraySort (peaks, pmPeaksCompareDescend);
    324324        clump = peaks->data[0];
    325         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
     325        psTrace (__func__, 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
    326326
    327327        // define section window for clump
    … …  
    369369        psfClump.dY = stats->clippedStdev;
    370370
    371         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
    372         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
    373         // these values should be pushed on the metadata somewhere
     371        psTrace (__func__, 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
     372        psTrace (__func__, 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
     373        // these values should be pushed on the recipe somewhere
    374374
    375375        psFree (stats);
    … …  
    379379    }
    380380
    381     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
     381    psTrace(__func__, 5, "---- end ----\n");
    382382    return (psfClump);
    383383}
    384384
    385385/******************************************************************************
    386     pmSourceRoughClass(source, metadata): make a guess at the source
     386    pmSourceRoughClass(source, recipe): make a guess at the source
    387387    classification.
    388  
    389     XXX: push the clump info into the metadata?
    390  
    391388    XXX: How can this function ever return FALSE?
    392  
    393     EAM: I moved S/N calculation to pmSourceMoments, using weight image
    394389*****************************************************************************/
    395390
    396 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
    397 {
    398     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
     391bool pmSourceRoughClass(psArray *sources, psMetadata *recipe, pmPSFClump clump)
     392{
     393    psTrace(__func__, 5, "---- begin ----");
    399394
    400395    psBool rc = true;
    … …  
    406401    int Ncr      = 0;
    407402    int Nsatstar = 0;
    408     // psRegion allArray = psRegionSet (0, 0, 0, 0);
    409403    psRegion inner;
    410404
    … …  
    412406    psVector *starsn = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    413407
    414     // check return status value (do these exist?)
     408    // get basic parameters, or set defaults
    415409    bool status;
    416     psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
    417     psF32 PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, metadata, "PSF_CLUMP_NSIGMA");
     410    float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
     411    if (!status)
     412        PSF_SN_LIM = 20.0;
     413    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    418414    if (!status)
    419415        PSF_CLUMP_NSIGMA = 1.5;
     416    float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    420417
    421418    // XXX allow clump size to be scaled relative to sigmas?
    … …  
    423420    for (psS32 i = 0 ; i < sources->n ; i++) {
    424421
    425         pmSource *tmpSrc = (pmSource *) sources->data[i];
    426 
    427         tmpSrc->peak->type = 0;
    428 
    429         psF32 sigX = tmpSrc->moments->Sx;
    430         psF32 sigY = tmpSrc->moments->Sy;
     422        pmSource *source = (pmSource *) sources->data[i];
     423
     424        source->peak->type = 0;
     425
     426        psF32 sigX = source->moments->Sx;
     427        psF32 sigY = source->moments->Sy;
    431428
    432429        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    433         // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
    434         inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2);
    435         inner = psRegionForImage (tmpSrc->mask, inner);
    436         int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PM_MASK_SAT);
     430        inner = psRegionForSquare (source->peak->x, source->peak->y, 2);
     431        inner = psRegionForImage (source->mask, inner);
     432        int Nsatpix = psImageCountPixelMask (source->mask, inner, PM_MASK_SAT);
    437433
    438434        // saturated star (size consistent with PSF or larger)
    … …  
    441437        big = true;
    442438        if ((Nsatpix > 1) && big) {
    443             tmpSrc->type = PM_SOURCE_TYPE_STAR;
    444             tmpSrc->mode |= PM_SOURCE_MODE_SATSTAR;
     439            source->type = PM_SOURCE_TYPE_STAR;
     440            source->mode |= PM_SOURCE_MODE_SATSTAR;
     441            // recalculate moments here with larger box?
     442            pmSourceMoments (source, INNER_RADIUS);
    445443            Nsatstar ++;
    446444            continue;
    … …  
    449447        // saturated object (not a star, eg bleed trails, hot pixels)
    450448        if (Nsatpix > 1) {
    451             tmpSrc->type = PM_SOURCE_TYPE_SATURATED;
    452             tmpSrc->mode |= PM_SOURCE_MODE_DEFAULT;
     449            source->type = PM_SOURCE_TYPE_SATURATED;
     450            source->mode |= PM_SOURCE_MODE_DEFAULT;
    453451            Nsat ++;
    454452            continue;
    … …  
    460458        // XXX these limits are quite arbitrary
    461459        if ((sigX < 0.05) || (sigY < 0.05)) {
    462             tmpSrc->type = PM_SOURCE_TYPE_DEFECT;
    463             tmpSrc->mode |= PM_SOURCE_MODE_DEFAULT;
     460            source->type = PM_SOURCE_TYPE_DEFECT;
     461            source->mode |= PM_SOURCE_MODE_DEFAULT;
    464462            Ncr ++;
    465463            continue;
    … …  
    468466        // likely unsaturated extended source (too large to be stellar)
    469467        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    470             tmpSrc->type = PM_SOURCE_TYPE_EXTENDED;
    471             tmpSrc->mode |= PM_SOURCE_MODE_DEFAULT;
     468            source->type = PM_SOURCE_TYPE_EXTENDED;
     469            source->mode |= PM_SOURCE_MODE_DEFAULT;
    472470            Next ++;
    473471            continue;
    … …  
    475473
    476474        // the rest are probable stellar objects
    477         starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
     475        starsn->data.F32[starsn->n] = source->moments->SN;
    478476        starsn->n ++;
    479477        Nstar ++;
    … …  
    481479        // PSF star (within 1.5 sigma of clump center, S/N > limit)
    482480        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    483         if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
    484             tmpSrc->type = PM_SOURCE_TYPE_STAR;
    485             tmpSrc->mode |= PM_SOURCE_MODE_PSFSTAR;
     481        if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
     482            source->type = PM_SOURCE_TYPE_STAR;
     483            source->mode |= PM_SOURCE_MODE_PSFSTAR;
    486484            Npsf ++;
    487485            continue;
    … …  
    489487
    490488        // random type of star
    491         tmpSrc->type = PM_SOURCE_TYPE_STAR;
    492         tmpSrc->mode |= PM_SOURCE_MODE_DEFAULT;
    493     }
    494 
    495     {
    496         psStats *stats  = NULL;
    497         stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    498         stats = psVectorStats (stats, starsn, NULL, NULL, 0);
    499         psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
    500         psFree (stats);
    501         psFree (starsn);
    502     }
    503 
    504     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
    505     psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
    506     psTrace (".pmObjects.pmSourceRoughClass", 2, "Next:     %3d\n", Next);
    507     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
    508     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
    509     psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);
    510 
    511     psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
     489        source->type = PM_SOURCE_TYPE_STAR;
     490        source->mode |= PM_SOURCE_MODE_DEFAULT;
     491    }
     492
     493    psStats *stats  = NULL;
     494    stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
     495    stats = psVectorStats (stats, starsn, NULL, NULL, 0);
     496    psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
     497    psFree (stats);
     498    psFree (starsn);
     499
     500    psTrace (__func__, 2, "Nstar:    %3d\n", Nstar);
     501    psTrace (__func__, 2, "Npsf:     %3d\n", Npsf);
     502    psTrace (__func__, 2, "Next:     %3d\n", Next);
     503    psTrace (__func__, 2, "Nsatstar: %3d\n", Nsatstar);
     504    psTrace (__func__, 2, "Nsat:     %3d\n", Nsat);
     505    psTrace (__func__, 2, "Ncr:      %3d\n", Ncr);
     506
     507    psTrace(__func__, 5, "---- end ----\n");
    512508    return(rc);
    513509}
    … …  
    535531                     psF32 radius)
    536532{
    537     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
     533    psTrace(__func__, 5, "---- begin ----\n");
    538534    PS_ASSERT_PTR_NON_NULL(source, false);
    539535    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    … …  
    641637    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    642638    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    643         psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
    644         psTrace("psModules.objects", 3, "---- %s(false) end ----\n", __func__);
     639        psTrace (__func__, 3, "no valid pixels for source\n");
     640        psTrace(__func__, 5, "---- end (false) ----\n");
    645641        return (false);
    646642    }
    647643
    648     psTrace (".psModules.pmSourceMoments", 5,
    649              "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
     644    psTrace (__func__, 4, "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
    650645             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
    651646
    … …  
    658653    y = Y1/Sum;
    659654    if ((fabs(x) > radius) || (fabs(y) > radius)) {
    660         psTrace (".psModules.pmSourceMoments", 3,
    661                  "large centroid swing; invalid peak %d, %d\n",
     655        psTrace (__func__, 3, "large centroid swing; invalid peak %d, %d\n",
    662656                 source->peak->x, source->peak->y);
    663         psTrace("psModules.objects", 3, "---- %s(false) end ----\n", __func__);
     657        psTrace(__func__, 5, "---- end(false)  ----\n");
    664658        return (false);
    665659    }
    … …  
    679673    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
    680674
    681     psTrace (".psModules.pmSourceMoments", 4,
     675    psTrace (__func__, 4,
    682676             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
    683677             sky, Sum, source->moments->x, source->moments->y,
    684678             source->moments->Sx, source->moments->Sy, source->moments->Sxy);
    685679
    686     psTrace("psModules.objects", 3, "---- %s(true) end ----\n", __func__);
     680    psTrace(__func__, 5, "---- end ----\n");
    687681    return(true);
    688682}
Note: See TracChangeset for help on using the changeset viewer.