IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 19, 2009, 7:59:50 AM (17 years ago)
Author:
beaumont
Message:

Added visualizations to ppSub. Set up a single variable in pmVisual to control when plots are drawn.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branch_20090215/psphot/src/psphotSourceStats.c

    r21359 r21536  
    11# include "psphotInternal.h"
    2 
    3 bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe);
    42
    53psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
    … …  
    8078        for (int j = 0; j < cells->n; j++) {
    8179
    82             if (nThreads) {
    83                 // allocate a job -- if threads are not defined, this just runs the job
    84                 psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
    85 
    86                 psArrayAdd(job->args, 1, cells->data[j]); // sources
    87                 psArrayAdd(job->args, 1, recipe);
    88                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
    89                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    90 
    91                 if (!psThreadJobAddPending(job)) {
    92                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    93                     psFree (job);
    94                     return NULL;
    95                 }
    96                 psFree(job);
    97             } else {
     80            // allocate a job -- if threads are not defined, this just runs the job
     81            psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
     82
     83            psArrayAdd(job->args, 1, cells->data[j]); // sources
     84            psArrayAdd(job->args, 1, recipe);
     85            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
     86            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     87
     88            if (!psThreadJobAddPending(job)) {
     89                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     90                psFree (job);
     91                return NULL;
     92            }
     93            psFree(job);
     94
     95# if (0)
    9896                int nfail = 0;
    9997                int nmoments = 0;
    … …  
    104102                Nfail += nfail;
    105103                Nmoments += nmoments;
     104# endif
     105        }
     106
     107        // wait for the threads to finish and manage results
     108        if (!psThreadPoolWait (false)) {
     109            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     110            return NULL;
     111        }
     112
     113        // we have only supplied one type of job, so we can assume the types here
     114        psThreadJob *job = NULL;
     115        while ((job = psThreadJobGetDone()) != NULL) {
     116            if (job->args->n < 1) {
     117                fprintf (stderr, "error with job\n");
     118            } else {
     119                psScalar *scalar = NULL;
     120                scalar = job->args->data[2];
     121                Nmoments += scalar->data.S32;
     122                scalar = job->args->data[3];
     123                Nfail += scalar->data.S32;
    106124            }
    107         }
    108 
    109         if (nThreads) {
    110             // wait for the threads to finish and manage results
    111             if (!psThreadPoolWait (false)) {
    112                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    113                 return NULL;
    114             }
    115 
    116             // we have only supplied one type of job, so we can assume the types here
    117             psThreadJob *job = NULL;
    118             while ((job = psThreadJobGetDone()) != NULL) {
    119                 if (job->args->n < 1) {
    120                     fprintf (stderr, "error with job\n");
    121                 } else {
    122                     psScalar *scalar = NULL;
    123                     scalar = job->args->data[2];
    124                     Nmoments += scalar->data.S32;
    125                     scalar = job->args->data[3];
    126                     Nfail += scalar->data.S32;
    127                 }
    128                 psFree(job);
    129             }
     125            psFree(job);
    130126        }
    131127    }
    … …  
    148144    psArray *sources                = job->args->data[0];
    149145    psMetadata *recipe              = job->args->data[1];
     146
     147    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     148    if (!status) return false;
     149    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
     150    if (!status) return false;
     151    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     152    if (!status) return false;
     153
     154    // bit-masks to test for good/bad pixels
     155    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
     156    assert (maskVal);
     157
     158    // bit-mask to mark pixels not used in analysis
     159    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     160    assert (markVal);
     161
     162    // maskVal is used to test for rejected pixels, and must include markVal
     163    maskVal |= markVal;
     164
     165    // threaded measurement of the sources moments
     166    int Nfail = 0;
     167    int Nmoments = 0;
     168    for (int i = 0; i < sources->n; i++) {
     169        pmSource *source = sources->data[i];
     170
     171        // skip faint sources for moments measurement
     172        if (source->peak->SN < MIN_SN) {
     173            source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN;
     174            continue;
     175        }
     176
     177        // measure a local sky value
     178        // the local sky is now ignored; kept here for reference only
     179        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
     180        if (!status) {
     181            source->mode |= PM_SOURCE_MODE_SKY_FAILURE;
     182            psErrorClear(); // XXX re-consider the errors raised here
     183            Nfail ++;
     184            continue;
     185        }
     186
     187        // measure the local sky variance (needed if noise is not sqrt(signal))
     188        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
     189        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
     190        if (!status) {
     191            source->mode |= PM_SOURCE_MODE_SKYVAR_FAILURE;
     192            Nfail ++;
     193            psErrorClear();
     194            continue;
     195        }
     196
     197        // measure basic source moments
     198        status = pmSourceMoments (source, RADIUS);
     199        if (status) {
     200            Nmoments ++;
     201            continue;
     202        }
     203
     204        // if no valid pixels, or massive swing, likely saturated source,
     205        // try a much larger box
     206        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
     207        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
     208        status = pmSourceMoments (source, BIG_RADIUS);
     209        if (status) {
     210            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
     211            Nmoments ++;
     212            continue;
     213        }
     214
     215        source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
     216        Nfail ++;
     217        psErrorClear();
     218        continue;
     219    }
     220
     221    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     222    scalar = job->args->data[2];
     223    scalar->data.S32 = Nmoments;
     224
     225    scalar = job->args->data[3];
     226    scalar->data.S32 = Nfail;
     227   
     228    return true;
     229}
     230
     231# if (0)
     232bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) {
     233
     234    bool status = false;
     235    float BIG_RADIUS;
    150236
    151237    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    … …  
    219305
    220306    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    221     scalar = job->args->data[2];
    222     scalar->data.S32 = Nmoments;
    223 
    224     scalar = job->args->data[3];
    225     scalar->data.S32 = Nfail;
    226    
    227     return true;
    228 }
    229 
    230 bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) {
    231 
    232     bool status = false;
    233     float BIG_RADIUS;
    234 
    235     float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    236     if (!status) return false;
    237     float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    238     if (!status) return false;
    239     float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    240     if (!status) return false;
    241 
    242     // bit-masks to test for good/bad pixels
    243     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    244     assert (maskVal);
    245 
    246     // bit-mask to mark pixels not used in analysis
    247     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    248     assert (markVal);
    249 
    250     // maskVal is used to test for rejected pixels, and must include markVal
    251     maskVal |= markVal;
    252 
    253     // threaded measurement of the sources moments
    254     int Nfail = 0;
    255     int Nmoments = 0;
    256     for (int i = 0; i < sources->n; i++) {
    257         pmSource *source = sources->data[i];
    258 
    259         // skip faint sources for moments measurement
    260         if (source->peak->SN < MIN_SN) {
    261             continue;
    262         }
    263 
    264         // measure a local sky value
    265         // the local sky is now ignored; kept here for reference only
    266         status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    267         if (!status) {
    268             psErrorClear(); // XXX re-consider the errors raised here
    269             Nfail ++;
    270             continue;
    271         }
    272 
    273         // measure the local sky variance (needed if noise is not sqrt(signal))
    274         // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
    275         status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    276         if (!status) {
    277             Nfail ++;
    278             psErrorClear();
    279             continue;
    280         }
    281 
    282         // measure basic source moments
    283         status = pmSourceMoments (source, RADIUS);
    284         if (status) {
    285             Nmoments ++;
    286             continue;
    287         }
    288 
    289         // if no valid pixels, or massive swing, likely saturated source,
    290         // try a much larger box
    291         BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    292         psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    293         status = pmSourceMoments (source, BIG_RADIUS);
    294         if (status) {
    295             Nmoments ++;
    296             continue;
    297         }
    298 
    299         Nfail ++;
    300         psErrorClear();
    301         continue;
    302     }
    303 
    304     // change the value of a scalar on the array (wrap this and put it in psArray.h)
    305307    *nmoments = Nmoments;
    306308    *nfail = Nfail;
    … …  
    308310    return true;
    309311}
     312# endif
Note: See TracChangeset for help on using the changeset viewer.