IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28911


Ignore:
Timestamp:
Aug 13, 2010, 12:20:26 PM (16 years ago)
Author:
eugene
Message:

thread the extended source analysis; update pmSourcePhotometry API

Location:
branches/eam_branches/ipp-20100621/psphot/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psphot/src/psphot.h

    r28821 r28911  
    106106bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule);
    107107bool            psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     108bool            psphotExtendedSourceFits_Threaded (psThreadJob *job);
    108109
    109110bool            psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule);

  • branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c

    r28782 r28911  
    3131// non-linear model fitting for extended sources
    3232bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     33
     34    bool status;
     35    int Next = 0;
     36    int Nconvolve = 0;
     37    int NconvolvePass = 0;
     38    int Nplain = 0;
     39    int NplainPass = 0;
     40
     41    psTimerStart ("psphot.extended");
     42
     43    // find the currently selected readout
     44    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     45    psAssert (file, "missing file?");
     46
     47    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     48    psAssert (readout, "missing readout?");
     49
     50    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     51    psAssert (detections, "missing detections?");
     52
     53    psArray *sources = detections->allSources;
     54    psAssert (sources, "missing sources?");
     55
     56    if (!sources->n) {
     57        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     58        return true;
     59    }
     60
     61    // determine the number of allowed threads
     62    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     63    if (!status) {
     64        nThreads = 0;
     65    }
     66
     67    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     68    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     69    assert (maskVal);
     70
     71    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
     72    assert (markVal);
     73
     74    // maskVal is used to test for rejected pixels, and must include markVal
     75    maskVal |= markVal;
     76
     77    // select the collection of desired models
     78    psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
     79    if (!status) {
     80        psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");
     81        return true;
     82    }
     83    if (models->list->n == 0) {
     84        psWarning ("extended source model fits requested but no models are specified\n");
     85        return true;
     86    }
     87
     88    // validate the model entries
     89    psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
     90    psMetadataItem *item = NULL;
     91    while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
     92
     93      if (item->type != PS_DATA_METADATA) {
     94        psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
     95        // XXX we could cull the bad entries or build a validated model folder
     96      }
     97
     98      psMetadata *model = (psMetadata *) item->data.md;
     99
     100      // check on the model type
     101      char *modelName = psMetadataLookupStr (&status, model, "MODEL");
     102      int modelType = pmModelClassGetType (modelName);
     103      if (modelType < 0) {
     104        psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
     105      }
     106      psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
     107
     108      // check on the SNLIM, set a float value
     109      char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
     110      if (!status) {
     111        psAbort("SNLIM not defined for extended source model %s\n", item->name);
     112      }
     113      float SNlim = atof (SNword);
     114      psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
     115
     116      // check on the PSF-Convolution status
     117      char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
     118      if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
     119        psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
     120      }
     121      bool convolved = !strcasecmp (convolvedWord, "true");
     122      psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
     123    }
     124    psFree (iter);
     125
     126    // option to limit analysis to a specific region
     127    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     128    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
     129    *AnalysisRegion = psRegionForImage(readout->image, psRegionFromString (region));
     130    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
     131
     132    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
     133    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
     134    assert (status);
     135
     136    // source analysis is done in S/N order (brightest first)
     137    sources = psArraySort (sources, pmSourceSortBySN);
     138
     139    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     140    int Cx = 1, Cy = 1;
     141    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     142
     143    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     144
     145    for (int i = 0; i < cellGroups->n; i++) {
     146
     147        psArray *cells = cellGroups->data[i];
     148
     149        for (int j = 0; j < cells->n; j++) {
     150
     151            // allocate a job -- if threads are not defined, this just runs the job
     152            psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_FIT");
     153
     154            psArrayAdd(job->args, 1, readout);
     155            psArrayAdd(job->args, 1, cells->data[j]); // sources
     156            psArrayAdd(job->args, 1, models);
     157            psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer
     158
     159            PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32);
     160            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
     161            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
     162
     163            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     164            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nconvolve
     165            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NconvolvePass
     166            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain
     167            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
     168
     169            if (!psThreadJobAddPending(job)) {
     170                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     171                psFree(AnalysisRegion);
     172                return false;
     173            }
     174        }
     175
     176        // wait for the threads to finish and manage results
     177        if (!psThreadPoolWait (false)) {
     178            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     179            psFree(AnalysisRegion);
     180            return false;
     181        }
     182
     183        // we have only supplied one type of job, so we can assume the types here
     184        psThreadJob *job = NULL;
     185        while ((job = psThreadJobGetDone()) != NULL) {
     186            if (job->args->n < 1) {
     187                fprintf (stderr, "error with job\n");
     188            } else {
     189                psScalar *scalar = NULL;
     190                scalar = job->args->data[7];
     191                Next += scalar->data.S32;
     192                scalar = job->args->data[8];
     193                Nconvolve += scalar->data.S32;
     194                scalar = job->args->data[9];
     195                NconvolvePass += scalar->data.S32;
     196                scalar = job->args->data[10];
     197                Nplain += scalar->data.S32;
     198                scalar = job->args->data[11];
     199                NplainPass += scalar->data.S32;
     200            }
     201            psFree(job);
     202            }
     203    }
     204    psFree (cellGroups);
     205    psFree(AnalysisRegion);
     206
     207    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     208    psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
     209    psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
     210    return true;
     211}
     212
     213// non-linear model fitting for extended sources
     214bool psphotExtendedSourceFits_Threaded (psThreadJob *job) {
    33215
    34216    bool status;
    … …  
    40222    bool savePics = false;
    41223    float radius;
    42 
    43     // find the currently selected readout
    44     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    45     psAssert (file, "missing file?");
    46 
    47     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    48     psAssert (readout, "missing readout?");
    49 
    50     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    51     psAssert (detections, "missing detections?");
    52 
    53     psArray *sources = detections->allSources;
    54     psAssert (sources, "missing sources?");
    55 
    56     if (!sources->n) {
    57         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
    58         return true;
    59     }
    60 
    61     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    62     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    63     assert (maskVal);
    64 
    65     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
    66     assert (markVal);
    67 
    68     // maskVal is used to test for rejected pixels, and must include markVal
    69     maskVal |= markVal;
    70 
    71     // select the collection of desired models
    72     psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
    73     if (!status) {
    74         psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");
    75         return true;
    76     }
    77     if (models->list->n == 0) {
    78         psWarning ("extended source model fits requested but no models are specified\n");
    79         return true;
    80     }
    81 
    82     // validate the model entries
    83     psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
    84     psMetadataItem *item = NULL;
    85     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
    86 
    87       if (item->type != PS_DATA_METADATA) {
    88         psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
    89         // XXX we could cull the bad entries or build a validated model folder
    90       }
    91 
    92       psMetadata *model = (psMetadata *) item->data.md;
    93 
    94       // check on the model type
    95       char *modelName = psMetadataLookupStr (&status, model, "MODEL");
    96       int modelType = pmModelClassGetType (modelName);
    97       if (modelType < 0) {
    98         psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
    99       }
    100       psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
    101 
    102       // check on the SNLIM, set a float value
    103       char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
    104       if (!status) {
    105         psAbort("SNLIM not defined for extended source model %s\n", item->name);
    106       }
    107       float SNlim = atof (SNword);
    108       psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
    109 
    110       // check on the PSF-Convolution status
    111       char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
    112       if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
    113         psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
    114       }
    115       bool convolved = !strcasecmp (convolvedWord, "true");
    116       psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
    117     }
    118     psFree (iter);
    119 
    120     // option to limit analysis to a specific region
    121     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    122     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    123     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    124 
    125     // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
    126     int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
    127     assert (status);
    128 
    129     // source analysis is done in S/N order (brightest first)
    130     sources = psArraySort (sources, pmSourceSortBySN);
     224    psScalar *scalar = NULL;
     225
     226    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
     227    pmReadout *readout      = job->args->data[0];
     228    psArray *sources        = job->args->data[1];
     229    psMetadata *models      = job->args->data[2];
     230    psRegion *region        = job->args->data[3];
     231    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     232    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     233    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
    131234
    132235    // Define source fitting parameters for extended source fits
    … …  
    150253
    151254        // XXX this should use peak?
    152         if (source->peak->x < AnalysisRegion.x0) continue;
    153         if (source->peak->y < AnalysisRegion.y0) continue;
    154         if (source->peak->x > AnalysisRegion.x1) continue;
    155         if (source->peak->y > AnalysisRegion.y1) continue;
     255        if (source->peak->x < region->x0) continue;
     256        if (source->peak->y < region->y0) continue;
     257        if (source->peak->x > region->x1) continue;
     258        if (source->peak->y > region->y1) continue;
    156259
    157260        // if model is NULL, we don't have a starting guess
    … …  
    166269
    167270        // set the radius based on the footprint (also sets the mask pixels)
    168         if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
     271        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
     272            psFree (fitOptions)
     273            return false;
     274        }
    169275
    170276        // XXX note that this changes the source moments that are published...
    … …  
    339445    psFree (fitOptions);
    340446
    341     psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
    342     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
    343     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
     447    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     448    scalar = job->args->data[7];
     449    scalar->data.S32 = Next;
     450
     451    scalar = job->args->data[8];
     452    scalar->data.S32 = Nconvolve;
     453
     454    scalar = job->args->data[9];
     455    scalar->data.S32 = NconvolvePass;
     456
     457    scalar = job->args->data[10];
     458    scalar->data.S32 = Nplain;
     459
     460    scalar = job->args->data[11];
     461    scalar->data.S32 = NplainPass;
     462
    344463    return true;
    345464}

  • branches/eam_branches/ipp-20100621/psphot/src/psphotLoadSRCTEXT.c

    r25983 r28911  
    8484            source->peak->yf   = PAR[PM_PAR_YPOS]; // but we know the pixel coordinate
    8585
    86             source->pixWeight = 1.0;
     86            source->pixWeightNotBad = 1.0;
     87            source->pixWeightNotPoor = 1.0;
    8788            source->crNsigma  = 0.0;
    8889            source->extNsigma = 0.0;

  • branches/eam_branches/ipp-20100621/psphot/src/psphotMagnitudes.c

    r28905 r28911  
    259259            psArrayAdd(job->args, 1, cells->data[j]); // sources
    260260            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     261            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
    261262
    262263            if (!psThreadJobAddPending(job)) {
    … …  
    295296    psArray *sources                = job->args->data[0];
    296297    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[1],PS_TYPE_IMAGE_MASK_DATA);
     298    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
    297299
    298300    for (int i = 0; i < sources->n; i++) {
    … …  
    303305        if (model == NULL) {
    304306          psTrace ("psphot", 3, "fail mag : no valid model");
    305           source->pixWeight = NAN;
     307          source->pixWeightNotBad = NAN;
     308          source->pixWeightNotPoor = NAN;
    306309          continue;
    307310        }
    308311
    309         status = pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
     312        status = pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal);
    310313        if (!status) {
    311314          psTrace ("psphot", 3, "fail to measure pixel weight");
    312           source->pixWeight = NAN;
     315          source->pixWeightNotBad = NAN;
     316          source->pixWeightNotPoor = NAN;
    313317          continue;
    314318        }

  • branches/eam_branches/ipp-20100621/psphot/src/psphotSetThreads.c

    r28643 r28911  
    1515    psFree(task);
    1616
    17     task = psThreadTaskAlloc("PSPHOT_PSF_WEIGHTS", 2);
     17    task = psThreadTaskAlloc("PSPHOT_PSF_WEIGHTS", 3);
    1818    task->function = &psphotPSFWeights_Threaded;
    1919    psThreadTaskAdd(task);
    … …  
    3535    psFree(task);
    3636
     37    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 12);
     38    task->function = &psphotExtendedSourceFits_Threaded;
     39    psThreadTaskAdd(task);
     40    psFree(task);
     41
    3742    return true;
    3843}

  • branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c

    r28702 r28911  
    370370                 source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY);
    371371
    372         // partially-masked sources are more likely to be mis-measured PSFs
     372        // XXX double check on ths stuff!! partially-masked sources are more likely to be mis-measured PSFs
    373373        float sizeBias = 1.0;
    374         if (source->pixWeight < 0.9) {
     374        if (source->pixWeightNotBad < 0.9) {
     375            sizeBias = 3.0;
     376        }
     377        if (source->pixWeightNotPoor < 0.9) {
    375378            sizeBias = 3.0;
    376379        }
    … …  
    407410        if (isCR) {
    408411            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
    409                     source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     412                    source->peak->xf,source->peak->yf,source->pixWeightNotBad,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    410413                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    411414            source->mode |= PM_SOURCE_MODE_DEFECT;
Note: See TracChangeset for help on using the changeset viewer.