IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 26, 2010, 9:18:39 AM (16 years ago)
Author:
Serge CHASTEL
Message:

Merging trunk in branch

Location:
branches/sc_branches/trunkTest
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/sc_branches/trunkTest

  • branches/sc_branches/trunkTest/psphot

    • Property svn:mergeinfo deleted

  • branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceFits.c

    r28013 r29060  
    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;
    … …  
    39221    int NplainPass = 0;
    40222    bool savePics = false;
    41 
    42     // find the currently selected readout
    43     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    44     psAssert (file, "missing file?");
    45 
    46     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    47     psAssert (readout, "missing readout?");
    48 
    49     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    50     psAssert (detections, "missing detections?");
    51 
    52     psArray *sources = detections->allSources;
    53     psAssert (sources, "missing sources?");
    54 
    55     if (!sources->n) {
    56         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
    57         return true;
    58     }
    59 
    60     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    61     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    62     assert (maskVal);
    63 
    64     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
    65     assert (markVal);
    66 
    67     // maskVal is used to test for rejected pixels, and must include markVal
    68     maskVal |= markVal;
    69 
    70     // select the collection of desired models
    71     psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
    72     if (!status) {
    73         psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");
    74         return true;
    75     }
    76     if (models->list->n == 0) {
    77         psWarning ("extended source model fits requested but no models are specified\n");
    78         return true;
    79     }
    80 
    81     // validate the model entries
    82     psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
    83     psMetadataItem *item = NULL;
    84     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
    85 
    86       if (item->type != PS_DATA_METADATA) {
    87         psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
    88         // XXX we could cull the bad entries or build a validated model folder
    89       }
    90 
    91       psMetadata *model = (psMetadata *) item->data.md;
    92 
    93       // check on the model type
    94       char *modelName = psMetadataLookupStr (&status, model, "MODEL");
    95       int modelType = pmModelClassGetType (modelName);
    96       if (modelType < 0) {
    97         psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
    98       }
    99       psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
    100 
    101       // check on the SNLIM, set a float value
    102       char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
    103       if (!status) {
    104         psAbort("SNLIM not defined for extended source model %s\n", item->name);
    105       }
    106       float SNlim = atof (SNword);
    107       psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
    108 
    109       // check on the PSF-Convolution status
    110       char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
    111       if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
    112         psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
    113       }
    114       bool convolved = !strcasecmp (convolvedWord, "true");
    115       psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
    116     }
    117     psFree (iter);
    118 
    119     // option to limit analysis to a specific region
    120     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    121     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    122     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    123 
    124     // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
    125     int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
    126     assert (status);
    127 
    128     // source analysis is done in S/N order (brightest first)
    129     sources = psArraySort (sources, pmSourceSortBySN);
     223    float radius;
     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);
     234
     235    // Define source fitting parameters for extended source fits
     236    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
     237    fitOptions->mode          = PM_SOURCE_FIT_EXT;
     238    // XXX for now, use the defaults for the rest:
     239    // fitOptions->nIter         = fitIter;
     240    // fitOptions->tol           = fitTol;
     241    // fitOptions->poissonErrors = poisson;
     242    // fitOptions->weight        = PS_SQR(skySig);
    130243
    131244    // choose the sources of interest
    … …  
    140253
    141254        // XXX this should use peak?
    142         if (source->peak->x < AnalysisRegion.x0) continue;
    143         if (source->peak->y < AnalysisRegion.y0) continue;
    144         if (source->peak->x > AnalysisRegion.x1) continue;
    145         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;
    146259
    147260        // if model is NULL, we don't have a starting guess
    … …  
    155268        Next ++;
    156269
     270        // set the radius based on the footprint (also sets the mask pixels)
     271        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
     272            fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
     273            psFree (fitOptions)
     274            return false;
     275        }
     276
     277        // XXX note that this changes the source moments that are published...
     278        // recalculate the source moments using the larger extended-source moments radius
     279        // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
     280        // this uses the footprint to judge both radius and aperture?
     281        // XXX save the psf-based moments for output
     282        if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
     283            fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
     284            // subtract the best fit from the object, leave local sky
     285            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     286            // XXX raise an error of some kind
     287            continue;
     288        }
     289
    157290        // save the modelFlux here in case we need to subtract it (for failure)
    158291        psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
    159 
     292        if (!modelFluxStart) {
     293            pmSourceCacheModel (source, maskVal);
     294            modelFluxStart = psMemIncrRefCounter (source->modelFlux);
     295            if (!modelFluxStart) {
     296                fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
     297                // XXX raise an error of some kind?
     298                continue;
     299            }
     300        }
     301       
    160302        if (savePics) {
    161303          psphotSaveImage (NULL, readout->image, "image.xp.fits");
    … …  
    183325          assert (status);
    184326
     327          // fprintf (stderr, "xfit: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN);
     328
    185329          // limit selection to some SN limit
    186330          assert (source->peak); // how can a source not have a peak?
    … …  
    201345          pmModel *modelFit = NULL;
    202346          if (convolved) {
    203               modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
     347              modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize);
    204348              if (!modelFit) {
    205349                  psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
    … …  
    215359              psFree (source->modelFlux);
    216360              source->modelFlux = NULL;
    217               modelFit = psphotFitEXT (readout, source, modelType, maskVal, markVal);
     361              modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal);
    218362              if (!modelFit) {
    219363                  psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
    … …  
    233377
    234378          // test for fit quality / result
     379          modelFit->fitRadius = radius;
    235380          psArrayAdd (source->modelFits, 4, modelFit);
    236381
    … …  
    311456        }
    312457    }
    313 
    314     psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
    315     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
    316     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
     458    psFree (fitOptions);
     459
     460    // fprintf (stderr, "xfit : tried %ld objects\n", sources->n);
     461
     462    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     463    scalar = job->args->data[7];
     464    scalar->data.S32 = Next;
     465
     466    scalar = job->args->data[8];
     467    scalar->data.S32 = Nconvolve;
     468
     469    scalar = job->args->data[9];
     470    scalar->data.S32 = NconvolvePass;
     471
     472    scalar = job->args->data[10];
     473    scalar->data.S32 = Nplain;
     474
     475    scalar = job->args->data[11];
     476    scalar->data.S32 = NplainPass;
     477
    317478    return true;
    318479}
Note: See TracChangeset for help on using the changeset viewer.