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/psphotBlendFit.c

    r21366 r21536  
    11# include "psphotInternal.h"
    2 
    3 bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources);
    42
    53// XXX I don't like this name
    … …  
    6260        for (int j = 0; j < cells->n; j++) {
    6361
    64             if (nThreads) {
    65                 // allocate a job -- if threads are not defined, this just runs the job
    66                 psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
    67                 psArray *newSources = psArrayAllocEmpty(16);
    68 
    69                 psArrayAdd(job->args, 1, readout);
    70                 psArrayAdd(job->args, 1, recipe);
    71                 psArrayAdd(job->args, 1, cells->data[j]); // sources
    72                 psArrayAdd(job->args, 1, psf);
    73                 psArrayAdd(job->args, 1, newSources); // return for new sources
    74                 psFree (newSources);
    75 
    76                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
    77                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
    78                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
    79                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    80 
    81                 if (!psThreadJobAddPending(job)) {
    82                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    83                     psFree (job);
    84                     return NULL;
    85                 }
    86                 psFree(job);
    87             } else {
     62            // allocate a job -- if threads are not defined, this just runs the job
     63            psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
     64            psArray *newSources = psArrayAllocEmpty(16);
     65
     66            psArrayAdd(job->args, 1, readout);
     67            psArrayAdd(job->args, 1, recipe);
     68            psArrayAdd(job->args, 1, cells->data[j]); // sources
     69            psArrayAdd(job->args, 1, psf);
     70            psArrayAdd(job->args, 1, newSources); // return for new sources
     71            psFree (newSources);
     72
     73            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
     74            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
     75            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     76            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     77
     78            if (!psThreadJobAddPending(job)) {
     79                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     80                psFree (job);
     81                return NULL;
     82            }
     83            psFree(job);
     84
     85# if (0)
     86            {
    8887                int nfit = 0;
    8988                int npsf = 0;
    … …  
    107106                psFree (newSources);
    108107            }
    109         }
    110 
    111         if (nThreads) {
    112             // wait for the threads to finish and manage results
    113             if (!psThreadPoolWait (false)) {
    114                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    115                 return NULL;
    116             }
    117 
    118             // we have only supplied one type of job, so we can assume the types here
    119             psThreadJob *job = NULL;
    120             while ((job = psThreadJobGetDone()) != NULL) {
    121                 if (job->args->n < 1) {
    122                     fprintf (stderr, "error with job\n");
    123                 } else {
    124                     psScalar *scalar = NULL;
    125                     scalar = job->args->data[5];
    126                     Nfit += scalar->data.S32;
    127                     scalar = job->args->data[6];
    128                     Npsf += scalar->data.S32;
    129                     scalar = job->args->data[7];
    130                     Next += scalar->data.S32;
    131                     scalar = job->args->data[8];
    132                     Nfail += scalar->data.S32;
    133 
    134                     // add these back onto sources
    135                     psArray *newSources = job->args->data[4];
    136                     for (int j = 0; j < newSources->n; j++) {
    137                         psArrayAdd (sources, 16, newSources->data[j]);
    138                     }
    139                 }
    140                 psFree(job);
    141             }
    142         }
     108# endif
     109        }
     110
     111        // wait for the threads to finish and manage results
     112        if (!psThreadPoolWait (false)) {
     113            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     114            return NULL;
     115        }
     116
     117        // we have only supplied one type of job, so we can assume the types here
     118        psThreadJob *job = NULL;
     119        while ((job = psThreadJobGetDone()) != NULL) {
     120            if (job->args->n < 1) {
     121                fprintf (stderr, "error with job\n");
     122            } else {
     123                psScalar *scalar = NULL;
     124                scalar = job->args->data[5];
     125                Nfit += scalar->data.S32;
     126                scalar = job->args->data[6];
     127                Npsf += scalar->data.S32;
     128                scalar = job->args->data[7];
     129                Next += scalar->data.S32;
     130                scalar = job->args->data[8];
     131                Nfail += scalar->data.S32;
     132
     133                // add these back onto sources
     134                psArray *newSources = job->args->data[4];
     135                for (int j = 0; j < newSources->n; j++) {
     136                    psArrayAdd (sources, 16, newSources->data[j]);
     137                }
     138            }
     139            psFree(job);
     140            }
    143141    }
    144142    psFree (cellGroups);
    … …  
    225223
    226224        // replace object in image
    227         if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
     225        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     226            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     227        }
     228        Nfit ++;
     229
     230        // try fitting PSFs or extended sources depending on source->mode
     231        // these functions subtract the resulting fitted source
     232        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     233            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
     234                source->type = PM_SOURCE_TYPE_EXTENDED;
     235                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
     236                Next ++;
     237                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     238                continue;
     239            }
     240        } else {
     241            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
     242                source->type = PM_SOURCE_TYPE_STAR;
     243                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
     244                Npsf ++;
     245                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     246                continue;
     247            }
     248        }
     249
     250        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
     251        Nfail ++;
     252
     253        // re-subtract the object, leave local sky
     254        pmSourceCacheModel (source, maskVal);
     255        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     256        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     257    }
     258
     259    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     260    scalar = job->args->data[5];
     261    scalar->data.S32 = Nfit;
     262
     263    scalar = job->args->data[6];
     264    scalar->data.S32 = Npsf;
     265
     266    scalar = job->args->data[7];
     267    scalar->data.S32 = Next;
     268
     269    scalar = job->args->data[8];
     270    scalar->data.S32 = Nfail;
     271
     272    return true;
     273}
     274
     275# if (0)
     276bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) {
     277
     278    bool status = false;
     279    int Nfit = 0;
     280    int Npsf = 0;
     281    int Next = 0;
     282    int Nfail = 0;
     283
     284    // bit-masks to test for good/bad pixels
     285    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
     286    assert (maskVal);
     287
     288    // bit-mask to mark pixels not used in analysis
     289    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     290    assert (markVal);
     291
     292    // maskVal is used to test for rejected pixels, and must include markVal
     293    maskVal |= markVal;
     294
     295    // S/N limit to perform full non-linear fits
     296    float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
     297
     298    // option to limit analysis to a specific region
     299    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     300    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     301    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
     302
     303    for (int i = 0; i < sources->n; i++) {
     304        pmSource *source = sources->data[i];
     305
     306        // skip non-astronomical objects (very likely defects)
     307        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
     308        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
     309        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     310        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     311
     312        // skip DBL second sources (ie, added by psphotFitBlob)
     313        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
     314
     315        // limit selection to some SN limit
     316        if (source->peak->SN < FIT_SN_LIM) continue;
     317
     318        // exclude sources outside optional analysis region
     319        if (source->peak->xf < AnalysisRegion.x0) continue;
     320        if (source->peak->yf < AnalysisRegion.y0) continue;
     321        if (source->peak->xf > AnalysisRegion.x1) continue;
     322        if (source->peak->yf > AnalysisRegion.y1) continue;
     323
     324        // if model is NULL, we don't have a starting guess
     325        if (source->modelPSF == NULL) continue;
     326
     327        // skip sources which are insignificant flux?
     328        // XXX this is somewhat ad-hoc
     329        if (source->modelPSF->params->data.F32[1] < 0.1) {
     330            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
     331                     source->modelPSF->params->data.F32[1],
     332                     source->modelPSF->params->data.F32[2],
     333                     source->modelPSF->params->data.F32[3]);
     334            continue;
     335        }
     336
     337        // replace object in image
     338        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    228339            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    229340        }
    … …  
    254365        pmSourceCacheModel (source, maskVal);
    255366        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    256         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
    257     }
    258 
    259     // change the value of a scalar on the array (wrap this and put it in psArray.h)
    260     scalar = job->args->data[5];
    261     scalar->data.S32 = Nfit;
    262 
    263     scalar = job->args->data[6];
    264     scalar->data.S32 = Npsf;
    265 
    266     scalar = job->args->data[7];
    267     scalar->data.S32 = Next;
    268 
    269     scalar = job->args->data[8];
    270     scalar->data.S32 = Nfail;
    271 
    272     return true;
    273 }
    274 
    275 bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) {
    276 
    277     bool status = false;
    278     int Nfit = 0;
    279     int Npsf = 0;
    280     int Next = 0;
    281     int Nfail = 0;
    282 
    283     // bit-masks to test for good/bad pixels
    284     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    285     assert (maskVal);
    286 
    287     // bit-mask to mark pixels not used in analysis
    288     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    289     assert (markVal);
    290 
    291     // maskVal is used to test for rejected pixels, and must include markVal
    292     maskVal |= markVal;
    293 
    294     // S/N limit to perform full non-linear fits
    295     float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
    296 
    297     // option to limit analysis to a specific region
    298     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    299     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    300     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    301 
    302     for (int i = 0; i < sources->n; i++) {
    303         pmSource *source = sources->data[i];
    304 
    305         // skip non-astronomical objects (very likely defects)
    306         if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
    307         if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
    308         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    309         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    310 
    311         // skip DBL second sources (ie, added by psphotFitBlob)
    312         if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
    313 
    314         // limit selection to some SN limit
    315         if (source->peak->SN < FIT_SN_LIM) continue;
    316 
    317         // exclude sources outside optional analysis region
    318         if (source->peak->xf < AnalysisRegion.x0) continue;
    319         if (source->peak->yf < AnalysisRegion.y0) continue;
    320         if (source->peak->xf > AnalysisRegion.x1) continue;
    321         if (source->peak->yf > AnalysisRegion.y1) continue;
    322 
    323         // if model is NULL, we don't have a starting guess
    324         if (source->modelPSF == NULL) continue;
    325 
    326         // skip sources which are insignificant flux?
    327         // XXX this is somewhat ad-hoc
    328         if (source->modelPSF->params->data.F32[1] < 0.1) {
    329             psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
    330                      source->modelPSF->params->data.F32[1],
    331                      source->modelPSF->params->data.F32[2],
    332                      source->modelPSF->params->data.F32[3]);
    333             continue;
    334         }
    335 
    336         // replace object in image
    337         if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
    338             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    339         }
    340         Nfit ++;
    341 
    342         // try fitting PSFs or extended sources depending on source->mode
    343         // these functions subtract the resulting fitted source
    344         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    345             if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
    346                 source->type = PM_SOURCE_TYPE_EXTENDED;
    347                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    348                 Next ++;
    349                 continue;
    350             }
    351         } else {
    352             if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
    353                 source->type = PM_SOURCE_TYPE_STAR;
    354                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    355                 Npsf ++;
    356                 continue;
    357             }
    358         }
    359 
    360         psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
    361         Nfail ++;
    362 
    363         // re-subtract the object, leave local sky
    364         pmSourceCacheModel (source, maskVal);
    365         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    366         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     367        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    367368    }
    368369
    … …  
    375376    return true;
    376377}
     378# endif
Note: See TracChangeset for help on using the changeset viewer.