IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 31, 2010, 5:01:05 PM (16 years ago)
Author:
eugene
Message:

updates relative to 20091201, fixes for all psphot variants

Location:
branches/eam_branches/psphot.stack.20100120
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/psphot.stack.20100120

  • branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c

    r26691 r26748  
    4242    psAssert (!detections->newSources, "new sources already defined?");
    4343
     44    // XXX TEST:
     45    if (detections->allSources) {
     46        psphotMaskCosmicRayFootprintCheck(detections->allSources);
     47    }
     48    if (detections->newSources) {
     49        psphotMaskCosmicRayFootprintCheck(detections->newSources);
     50    }
     51
    4452    // determine the number of allowed threads
    4553    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    … …  
    5765    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    5866    psAssert (status, "missing BREAK_POINT?");
     67
     68    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); psAssert (status, "missing SKY_INNER_RADIUS");
     69    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
     70
     71    // bit-masks to test for good/bad pixels
     72    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
     73    psAssert (maskVal, "missing MASK.PSPHOT");
     74
     75    // bit-mask to mark pixels not used in analysis
     76    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     77    psAssert (markVal, "missing MARK.PSPHOT");
    5978
    6079    psArray *peaks = detections->peaks;
    … …  
    112131    }
    113132
     133    // if we have measured the window, we will be saving the modified version of these recipe values on readout->analysis
     134    float SIGMA = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
     135    if (!status) {
     136        SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
     137    }
     138    float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
     139    if (!status) {
     140        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     141    }
     142
    114143    // threaded measurement of the source magnitudes
    115144    int Nfail = 0;
    … …  
    133162
    134163            psArrayAdd(job->args, 1, cells->data[j]); // sources
    135             psArrayAdd(job->args, 1, recipe);
     164
     165            PS_ARRAY_ADD_SCALAR(job->args, INNER,   PS_TYPE_F32);
     166            PS_ARRAY_ADD_SCALAR(job->args, MIN_SN,  PS_TYPE_F32);
     167            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
     168            PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
     169
     170            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
     171            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
     172
    136173            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
    137174            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    … …  
    152189            psFree(sources);
    153190            return false;
     191        }
     192
     193        // we have only supplied one type of job, so we can assume the types here
     194        psThreadJob *job = NULL;
     195        while ((job = psThreadJobGetDone()) != NULL) {
     196            if (job->args->n < 1) {
     197                fprintf (stderr, "error with job\n");
     198            } else {
     199                psScalar *scalar = NULL;
     200                scalar = job->args->data[7];
     201                Nmoments += scalar->data.S32;
     202                scalar = job->args->data[8];
     203                Nfail += scalar->data.S32;
     204                scalar = job->args->data[9];
     205                Nfaint += scalar->data.S32;
     206            }
     207            psFree(job);
     208        }
     209    }
     210
     211    psFree (cellGroups);
     212
     213    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
     214
     215    psphotVisualShowMoments (sources);
     216
     217    // save the new sources on the detection structure:
     218    detections->newSources = sources;
     219
     220
     221    if (detections->allSources) {
     222        psphotMaskCosmicRayFootprintCheck(detections->allSources);
     223    }
     224    if (detections->newSources) {
     225        psphotMaskCosmicRayFootprintCheck(detections->newSources);
     226    }
     227
     228    return true;
     229}
     230
     231// this function is currently only called by psphotCheckExtSources
     232bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) {
     233
     234    bool status = false;
     235
     236    psTimerStart ("psphot.stats");
     237
     238    // select the appropriate recipe information
     239    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     240    assert (recipe);
     241
     242    // determine the number of allowed threads
     243    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     244    if (!status) {
     245        nThreads = 0;
     246    }
     247
     248    // determine properties (sky, moments) of initial sources
     249    float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     250    if (!status) return false;
     251
     252    OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
     253
     254    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
     255    if (!status) return false;
     256
     257    for (int i = 0; i < sources->n; i++) {
     258
     259        pmSource *source = sources->data[i];
     260        if (!source->peak) continue; // XXX how can we have a peak-less source?
     261
     262        // allocate space for moments
     263        if (!source->moments) {
     264            source->moments = pmMomentsAlloc();
     265        }
     266
     267        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
     268        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
     269        source->peak->assigned = true;
     270    }
     271
     272    if (!strcasecmp (breakPt, "PEAKS")) {
     273        psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
     274        psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
     275        psphotVisualShowMoments (sources);
     276        return sources;
     277    }
     278
     279    // XXX how else could we get the window size in?
     280    if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
     281        psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
     282        return NULL;
     283    }
     284
     285    // threaded measurement of the source magnitudes
     286    int Nfail = 0;
     287    int Nmoments = 0;
     288    int Nfaint = 0;
     289
     290    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     291    int Cx = 1, Cy = 1;
     292    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     293
     294    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     295
     296    for (int i = 0; i < cellGroups->n; i++) {
     297
     298        psArray *cells = cellGroups->data[i];
     299
     300        for (int j = 0; j < cells->n; j++) {
     301
     302            // allocate a job -- if threads are not defined, this just runs the job
     303            psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
     304
     305            // XXX: this must match the above
     306            psArrayAdd(job->args, 1, cells->data[j]); // sources
     307            psArrayAdd(job->args, 1, recipe);
     308            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
     309            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     310            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
     311
     312            if (!psThreadJobAddPending(job)) {
     313                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
     314                psFree (job);
     315                return NULL;
     316            }
     317            psFree(job);
     318        }
     319
     320        // wait for the threads to finish and manage results
     321        if (!psThreadPoolWait (false)) {
     322            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
     323            return NULL;
    154324        }
    155325
    … …  
    178348    psphotVisualShowMoments (sources);
    179349
    180     // save the new sources on the detection structure:
    181     detections->newSources = sources;
    182     return true;
    183 }
    184 
    185 // XXX fix the return status of this function...
    186 bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) {
    187 
    188     bool status = false;
    189 
    190     psTimerStart ("psphot.stats");
    191 
    192     // select the appropriate recipe information
    193     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    194     assert (recipe);
    195 
    196     // determine the number of allowed threads
    197     int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    198     if (!status) {
    199         nThreads = 0;
    200     }
    201 
    202     // determine properties (sky, moments) of initial sources
    203     float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    204     if (!status) return false;
    205 
    206     OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
    207 
    208     char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    209     if (!status) return false;
    210 
    211     for (int i = 0; i < sources->n; i++) {
    212 
    213         pmSource *source = sources->data[i];
    214         if (!source->peak) continue; // XXX how can we have a peak-less source?
    215 
    216         // allocate space for moments
    217         if (!source->moments) {
    218             source->moments = pmMomentsAlloc();
    219         }
    220 
    221         // allocate image, weight, mask arrays for each peak (square of radius OUTER)
    222         pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
    223         source->peak->assigned = true;
    224     }
    225 
    226     if (!strcasecmp (breakPt, "PEAKS")) {
    227         psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
    228         psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
    229         psphotVisualShowMoments (sources);
    230         return sources;
    231     }
    232 
    233     // XXX how else could we get the window size in?
    234     if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
    235         psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
    236         return NULL;
    237     }
    238 
    239     // threaded measurement of the source magnitudes
    240     int Nfail = 0;
    241     int Nmoments = 0;
    242     int Nfaint = 0;
    243 
    244     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    245     int Cx = 1, Cy = 1;
    246     psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
    247 
    248     psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
    249 
    250     for (int i = 0; i < cellGroups->n; i++) {
    251 
    252         psArray *cells = cellGroups->data[i];
    253 
    254         for (int j = 0; j < cells->n; j++) {
    255 
    256             // allocate a job -- if threads are not defined, this just runs the job
    257             psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
    258 
    259             psArrayAdd(job->args, 1, cells->data[j]); // sources
    260             psArrayAdd(job->args, 1, recipe);
    261             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
    262             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    263             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
    264 
    265             if (!psThreadJobAddPending(job)) {
    266                 psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
    267                 psFree (job);
    268                 return NULL;
    269             }
    270             psFree(job);
    271         }
    272 
    273         // wait for the threads to finish and manage results
    274         if (!psThreadPoolWait (false)) {
    275             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
    276             return NULL;
    277         }
    278 
    279         // we have only supplied one type of job, so we can assume the types here
    280         psThreadJob *job = NULL;
    281         while ((job = psThreadJobGetDone()) != NULL) {
    282             if (job->args->n < 1) {
    283                 fprintf (stderr, "error with job\n");
    284             } else {
    285                 psScalar *scalar = NULL;
    286                 scalar = job->args->data[2];
    287                 Nmoments += scalar->data.S32;
    288                 scalar = job->args->data[3];
    289                 Nfail += scalar->data.S32;
    290                 scalar = job->args->data[4];
    291                 Nfaint += scalar->data.S32;
    292             }
    293             psFree(job);
    294         }
    295     }
    296 
    297     psFree (cellGroups);
    298 
    299     psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
    300 
    301     psphotVisualShowMoments (sources);
    302 
    303350    return (sources);
    304351}
    … …  
    311358
    312359    psArray *sources                = job->args->data[0];
    313     psMetadata *recipe              = job->args->data[1];
    314 
    315     float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    316     if (!status) return false;
    317     float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    318     if (!status) return false;
    319     float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    320     if (!status) return false;
    321     float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
    322     if (!status) return false;
    323 
    324     // bit-masks to test for good/bad pixels
    325     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    326     assert (maskVal);
    327 
    328     // bit-mask to mark pixels not used in analysis
    329     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    330     assert (markVal);
     360
     361    float INNER                     = PS_SCALAR_VALUE(job->args->data[1],F32);
     362    float MIN_SN                    = PS_SCALAR_VALUE(job->args->data[2],F32);
     363    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[3],F32);
     364    float SIGMA                     = PS_SCALAR_VALUE(job->args->data[4],F32);
     365
     366    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     367    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
    331368
    332369    // maskVal is used to test for rejected pixels, and must include markVal
    … …  
    392429
    393430    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    394     scalar = job->args->data[2];
     431    scalar = job->args->data[7];
    395432    scalar->data.S32 = Nmoments;
    396433
    397     scalar = job->args->data[3];
     434    scalar = job->args->data[8];
    398435    scalar->data.S32 = Nfail;
    399436
    400     scalar = job->args->data[4];
     437    scalar = job->args->data[9];
    401438    scalar->data.S32 = Nfaint;
    402439
    … …  
    410447    bool status;
    411448
    412     float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    413     if (!status) return false;
     449    float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
     450    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
     451    psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
    414452
    415453    // XXX move this to a config file?
    … …  
    441479
    442480        // choose a grid scale that is a fixed fraction of the psf sigma^2
    443         psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
    444         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
    445         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
     481        float PSF_CLUMP_GRID_SCALE = 0.1*PS_SQR(sigma[i]);
     482        float MOMENTS_SX_MAX = 2.0*PS_SQR(sigma[i]);
     483        float MOMENTS_SY_MAX = 2.0*PS_SQR(sigma[i]);
    446484
    447485        // determine the PSF parameters from the source moment values
    448         pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
     486        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
    449487        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
    450488
    … …  
    485523
    486524    // choose a grid scale that is a fixed fraction of the psf sigma^2
    487     psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
    488     psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
    489     psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
    490     psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
    491     psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
     525    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
     526    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
     527    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
     528    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
     529    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
    492530
    493531    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
Note: See TracChangeset for help on using the changeset viewer.