IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21519


Ignore:
Timestamp:
Feb 16, 2009, 12:34:24 PM (17 years ago)
Author:
eugene
Message:

added additional flag values, threaded functions always use the threaded version

Location:
trunk/psphot/src
Files:
1 deleted
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotAddNoise.c

    r21183 r21519  
    5050
    5151        // skip sources which were not subtracted
    52         if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
     52        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
    5353
    5454        // select appropriate model

  • trunk/psphot/src/psphotApResid.c

    r21359 r21519  
    11# include "psphotInternal.h"
    2 
    3 bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal);
    42
    53# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
    … …  
    9593        for (int j = 0; j < cells->n; j++) {
    9694
    97             if (nThreads) {
    98                 // allocate a job -- if threads are not defined, this just runs the job
    99                 psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS");
    100 
    101                 psArrayAdd(job->args, 1, cells->data[j]); // sources
    102                 psArrayAdd(job->args, 1, psf);
    103                 PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
    104                 PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    105 
    106                 PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
    107                 PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
    108 
    109                 if (!psThreadJobAddPending(job)) {
    110                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    111                     psFree (job);
    112                     return false;
    113                 }
    114                 psFree(job);
    115             } else {
     95            // allocate a job -- if threads are not defined, this just runs the job
     96            psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS");
     97
     98            psArrayAdd(job->args, 1, cells->data[j]); // sources
     99            psArrayAdd(job->args, 1, psf);
     100            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
     101            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     102
     103            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
     104            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
     105
     106            if (!psThreadJobAddPending(job)) {
     107                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     108                psFree (job);
     109                return false;
     110            }
     111            psFree(job);
     112
     113# if (0)
    116114                int nskip = 0;
    117115                int nfail = 0;
    … …  
    123121                Nskip += nskip;
    124122                Nfail += nfail;
    125             }           
     123# endif
    126124
    127125        }
    128126
    129         if (nThreads) {
    130             // wait for the threads to finish and manage results
    131             if (!psThreadPoolWait (false)) {
    132                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    133                 return false;
     127        // wait for the threads to finish and manage results
     128        if (!psThreadPoolWait (false)) {
     129            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     130            return false;
     131        }
     132
     133        // we have only supplied one type of job, so we can assume the types here
     134        psThreadJob *job = NULL;
     135        while ((job = psThreadJobGetDone()) != NULL) {
     136            if (job->args->n < 1) {
     137                fprintf (stderr, "error with job\n");
     138            } else {
     139                psScalar *scalar = NULL;
     140                scalar = job->args->data[4];
     141                Nskip += scalar->data.S32;
     142                scalar = job->args->data[5];
     143                Nfail += scalar->data.S32;
    134144            }
    135 
    136             // we have only supplied one type of job, so we can assume the types here
    137             psThreadJob *job = NULL;
    138             while ((job = psThreadJobGetDone()) != NULL) {
    139                 if (job->args->n < 1) {
    140                     fprintf (stderr, "error with job\n");
    141                 } else {
    142                     psScalar *scalar = NULL;
    143                     scalar = job->args->data[4];
    144                     Nskip += scalar->data.S32;
    145                     scalar = job->args->data[5];
    146                     Nfail += scalar->data.S32;
    147                 }
    148                 psFree(job);
    149             }
     145            psFree(job);
    150146        }
    151147    }
    … …  
    497493            continue;
    498494        }
     495        source->mode |= PM_SOURCE_MODE_AP_MAGS;
    499496    }
    500497
    … …  
    509506}
    510507
     508# if (0)
    511509bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
    512510
    … …  
    542540    return true;
    543541}
     542# endif

  • trunk/psphot/src/psphotBlendFit.c

    r21366 r21519  
    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

  • trunk/psphot/src/psphotDeblendSatstars.c

    r20453 r21519  
    111111       
    112112        // any peaks within this contour should be dropped (mark as blend, drop after loop is done)
    113         // also drop any peaks which are too close to this peal
     113        // also drop any peaks which are too close to this peak
    114114        psVector *xv = contour->data[0];
    115115        psVector *yv = contour->data[1];

  • trunk/psphot/src/psphotDetectReadout.c

    r21392 r21519  
    2222    // Generate the mask and weight images, including the user-defined analysis region of interest
    2323    psphotSetMaskAndWeight (config, readout, recipe);
    24     if (!strcasecmp (breakPt, "NOTHING")) {
    25         return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);
    26     }
    2724
    2825    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    … …  
    3835
    3936    // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF)
    40     pmDetections *detections = psphotDetectionsFromSources (config, inSources);
    41     if (!detections || !detections->peaks) {
    42         psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars");
    43         return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL);
    44     }
    45 
    46     // construct sources and measure basic stats
    47     psArray *sources = psphotSourceStats (config, readout, detections);
    48     if (!sources) return false;
    49     if (!strcasecmp (breakPt, "PEAKS")) {
    50         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    51     }
    52 
    53     // classify sources based on moments, brightness
    54     if (!psphotRoughClass (readout, sources, recipe, havePSF)) {
    55         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    56         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    57     }
     37    psphotSetSourceParams (config, sources, psf);
    5838
    5939    // calculate source magnitudes
    … …  
    6343    return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    6444}
     45
     46// deep in pmSourcePixelWeight, we need the following items for each location of interest:
     47// pmModel *model (should be a realized version of the PSF, can have unity normalization)
     48// psImage *mask  (local subimage of mask for source)
     49// psImageMaskType maskVal (from recipe & mask header)
     50
     51{
     52    pmModel *model = pmModelFromPSFforXY (psf, x, y, 1.0);
     53}
     54

  • trunk/psphot/src/psphotExtendedSourceAnalysis.c

    r21183 r21519  
    6161
    6262        // replace object in image
    63         if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
     63        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    6464            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    6565        }
    … …  
    7373                psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    7474                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    75                 source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     75                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    7676                continue;
    7777            }
     78            source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
    7879        }
    7980
    … …  
    8586                psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    8687                Nisophot ++;
     88                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    8789            }
    8890        }
    … …  
    9597                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    9698                Npetro ++;
     99                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    97100            }
    98101        }
    … …  
    105108                psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    106109                Nkron ++;
     110                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    107111            }
    108112        }
    … …  
    116120            psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    117121            Nannuli ++;
     122            source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    118123        }
    119124
    120125        // re-subtract the object, leave local sky
    121126        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    122         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     127        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    123128    }
    124129

  • trunk/psphot/src/psphotExtendedSourceFits.c

    r21183 r21519  
    110110
    111111        // replace object in image
    112         if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
     112        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    113113            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    114114        }
    … …  
    170170              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
    171171                  NconvolvePass ++;
     172                  source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
    172173              }
    173174          } else {
    … …  
    184185              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
    185186                  NplainPass ++;
     187                  source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
    186188              }
    187189          }
    … …  
    224226
    225227          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    226           source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     228          source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    227229
    228230          psFree (modelFluxes);
    … …  
    251253        // subtract the best fit from the object, leave local sky
    252254        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    253         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     255        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    254256
    255257        // the initial model flux is no longer needed

  • trunk/psphot/src/psphotFindPeaks.c

    r20453 r21519  
    2121    }
    2222
    23     // correct the peak values to S/N = sqrt(significance)
    24     // get the peak flux from the unsmoothed image
    25     // the peak pixel coords are guaranteed to be on the image
     23    // Convert the peak values to S/N = sqrt(significance).
     24    // Get the peak flux from the unsmoothed image.
     25    // Rescale the peak position errors using the peak variance
     26    // The peak pixel coords are guaranteed to be on the image
    2627    int row0 = readout->image->row0;
    2728    int col0 = readout->image->col0;
    … …  
    3031        peak->SN = sqrt(peak->value);
    3132        peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
     33        if (readout->variance && isfinite (peak->dx)) {
     34            peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
     35        }
     36        if (readout->variance && isfinite (peak->dy)) {
     37            peak->dy *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]);
     38        }
    3239    }
    3340

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r21366 r21519  
    6363        pmSource *source = sources->data[i];
    6464
     65        // turn this bit off and turn it on again if we pass this test
     66        source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
     67
    6568        // skip non-astronomical objects (very likely defects)
    6669        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    6770        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    6871
    69         // if (source->type == PM_SOURCE_TYPE_STAR &&
     72        // do not include CRs in the full ensemble fit
    7073        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
    7174
    7275        if (final) {
    73             if (source->mode &  PM_SOURCE_MODE_SUBTRACTED) continue;
     76            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
    7477        } else {
    75             if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
     78            if (source->mode & PM_SOURCE_MODE_BLEND) continue;
    7679        }
    7780
    … …  
    9194        if (y > AnalysisRegion.y1) continue;
    9295
     96        source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
    9397        psArrayAdd (fitSources, 100, source);
    9498    }
    … …  
    185189    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    186190
     191    // XXXX **** philosophical question:
     192    // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit:
     193    // should retain the chisq and errors from the intermediate non-linear fit?
     194    // the non-linear fit provides better values for the position errors, and for
     195    // extended sources, the shape errors
     196
    187197    // adjust I0 for fitSources and subtract
    188198    for (int i = 0; i < fitSources->n; i++) {
    … …  
    194204            psAbort("linear fitted source is nan");
    195205        }
     206
    196207        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
    197208        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
    … …  
    200211        // subtract object
    201212        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    202         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     213        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    203214    }
    204215    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    … …  
    207218    for (int i = 0; final && (i < fitSources->n); i++) {
    208219        pmSource *source = fitSources->data[i];
     220        if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
    209221        pmModel *model = pmSourceGetModel (NULL, source);
    210222        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);

  • trunk/psphot/src/psphotGuessModels.c

    r21381 r21519  
    2222}
    2323
    24 bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal);
    25 
    2624// construct an initial PSF model for each object
    2725bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
    … …  
    6765        for (int j = 0; j < cells->n; j++) {
    6866
    69             if (nThreads) {
    70                 // allocate a job -- if threads are not defined, this just runs the job
    71                 psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
    72                 psArrayAdd(job->args, 1, readout);
    73                 psArrayAdd(job->args, 1, cells->data[j]); // sources
    74                 psArrayAdd(job->args, 1, psf);
    75 
    76                 // XXX change these to use abstract mask type info
    77                 PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    78                 PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
    79 
    80                 if (!psThreadJobAddPending(job)) {
    81                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    82                     psFree (job);
    83                     return false;
    84                 }
    85                 psFree(job);
    86             } else {
     67            // allocate a job -- if threads are not defined, this just runs the job
     68            psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
     69            psArrayAdd(job->args, 1, readout);
     70            psArrayAdd(job->args, 1, cells->data[j]); // sources
     71            psArrayAdd(job->args, 1, psf);
     72
     73            // XXX change these to use abstract mask type info
     74            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     75            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
     76
     77            if (!psThreadJobAddPending(job)) {
     78                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     79                psFree (job);
     80                return false;
     81            }
     82            psFree(job);
     83
     84# if (0)               
    8785                if (!psphotGuessModel_Unthreaded (readout, cells->data[j], psf, maskVal, markVal)) {
    8886                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    8987                    return false;
    9088                }
    91             }
    92         }
    93 
    94         if (nThreads) {
    95             // wait for the threads to finish and manage results
    96             // wait here for the threaded jobs to finish
    97             // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
    98             if (!psThreadPoolWait (false)) {
    99                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    100                 return false;
    101             }
    102 
    103             // we have only supplied one type of job, so we can assume the types here
    104             psThreadJob *job = NULL;
    105             while ((job = psThreadJobGetDone()) != NULL) {
    106                 // we have no returned data from this operation
    107                 if (job->args->n < 1) {
    108                     fprintf (stderr, "error with job\n");
    109                 }
    110                 psFree(job);
    111             }
     89# endif
     90        }
     91
     92        // wait for the threads to finish and manage results
     93        // wait here for the threaded jobs to finish
     94        // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
     95        if (!psThreadPoolWait (false)) {
     96            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     97            return false;
     98        }
     99
     100        // we have only supplied one type of job, so we can assume the types here
     101        psThreadJob *job = NULL;
     102        while ((job = psThreadJobGetDone()) != NULL) {
     103            // we have no returned data from this operation
     104            if (job->args->n < 1) {
     105                fprintf (stderr, "error with job\n");
     106            }
     107            psFree(job);
    112108        }
    113109    }
    … …  
    116112    int nMiss = 0;
    117113    for (int i = 0; i < sources->n; i++) {
    118 
    119114        pmSource *source = sources->data[i];
    120         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    121             source->mode &= ~PM_SOURCE_MODE_EXT_LIMIT;
     115        if (source->tmpFlags & PM_SOURCE_TMPF_MODEL_GUESS) {
    122116            continue;
    123117        }
    124 
    125118        nMiss ++;
    126119    }
    127     psLogMsg ("psphot.models", 4, "failed to build models for %d objects\n", nMiss);
     120    psAssert (nMiss == 0, "failed to attempt to build models for %d objects\n", nMiss);
    128121
    129122    psFree (cellGroups);
    … …  
    148141        pmSource *source = sources->data[i];
    149142
    150         // XXXX this is just for a test: use this to mark sources for which the model is measured
    151         // check later that all are used.
    152         source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     143        // this is used to mark sources for which the model is measured. We check later that
     144        // all are used.
     145        source->tmpFlags |= PM_SOURCE_TMPF_MODEL_GUESS;
    153146
    154147        // skip non-astronomical objects (very likely defects)
    … …  
    185178        pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC X5
    186179        if (modelPSF == NULL) {
    187             psError(PSPHOT_ERR_PSF, false,
    188                     "Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
     180            psWarning ("Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
    189181                    source->peak->y, source->peak->x);
    190             //
    191             // Try the centre of the image
    192             //
     182
     183            // Try the center of the image
    193184            modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
    194185            modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
    195186            modelPSF = pmModelFromPSF (modelEXT, psf);
    196187            if (modelPSF == NULL) {
    197                 psError(PSPHOT_ERR_PSF, false,
    198                         "Failed to determine PSF model at centre of image");
     188                psError(PSPHOT_ERR_PSF, false, "Failed to determine PSF model at center of image");
    199189                psFree(modelEXT);
    200190                return false;
    201191            }
    202 
    203192            source->mode |= PM_SOURCE_MODE_BADPSF;
    204193        }
    … …  
    221210}
    222211
     212# if (0)
    223213// construct models only for sources in the specified region
    224214bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
    … …  
    301291    return true;
    302292}
     293# endif

  • trunk/psphot/src/psphotMagnitudes.c

    r21392 r21519  
    11# include "psphotInternal.h"
    2 
    3 bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal);
    42
    53bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
    … …  
    5957        for (int j = 0; j < cells->n; j++) {
    6058
    61             if (nThreads) {
    62                 // allocate a job -- if threads are not defined, this just runs the job
    63                 psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES");
    64 
    65                 psArrayAdd(job->args, 1, cells->data[j]); // sources
    66                 psArrayAdd(job->args, 1, psf);
    67                 psArrayAdd(job->args, 1, binning);
    68                 psArrayAdd(job->args, 1, backModel);
    69                 psArrayAdd(job->args, 1, backStdev);
    70 
    71                 PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
    72                 PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    73                 PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
    74 
    75                 if (!psThreadJobAddPending(job)) {
    76                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    77                     psFree (job);
    78                     return false;
    79                 }
    80                 psFree(job);
    81             } else {
     59            // allocate a job -- if threads are not defined, this just runs the job
     60            psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES");
     61
     62            psArrayAdd(job->args, 1, cells->data[j]); // sources
     63            psArrayAdd(job->args, 1, psf);
     64            psArrayAdd(job->args, 1, binning);
     65            psArrayAdd(job->args, 1, backModel);
     66            psArrayAdd(job->args, 1, backStdev);
     67
     68            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
     69            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     70            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
     71
     72            if (!psThreadJobAddPending(job)) {
     73                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     74                psFree (job);
     75                return false;
     76            }
     77            psFree(job);
     78
     79# if (0)
    8280                int nap = 0;
    8381                if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) {
    … …  
    8684                }
    8785                Nap += nap;
    88             }
    89         }
    90 
    91         if (nThreads) {
    92             // wait for the threads to finish and manage results
    93             if (!psThreadPoolWait (false)) {
    94                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    95                 return false;
    96             }
    97 
    98             // we have only supplied one type of job, so we can assume the types here
    99             psThreadJob *job = NULL;
    100             while ((job = psThreadJobGetDone()) != NULL) {
    101                 if (job->args->n < 1) {
    102                     fprintf (stderr, "error with job\n");
    103                 } else {
    104                     psScalar *scalar = job->args->data[7];
    105                     Nap += scalar->data.S32;
    106                 }
    107                 psFree(job);
    108             }
     86# endif
     87        }
     88
     89        // wait for the threads to finish and manage results
     90        if (!psThreadPoolWait (false)) {
     91            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     92            return false;
     93        }
     94
     95        // we have only supplied one type of job, so we can assume the types here
     96        psThreadJob *job = NULL;
     97        while ((job = psThreadJobGetDone()) != NULL) {
     98            if (job->args->n < 1) {
     99                fprintf (stderr, "error with job\n");
     100            } else {
     101                psScalar *scalar = job->args->data[7];
     102                Nap += scalar->data.S32;
     103            }
     104            psFree(job);
    109105        }
    110106    }
    … …  
    162158}
    163159
     160# if (0)
    164161bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
    165162
    … …  
    198195    return true;
    199196}
     197# endif
    200198
    201199bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources) {
    … …  
    295293        }
    296294
    297         status = pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal);
     295        status = pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
    298296        if (!status) {
    299297          psTrace ("psphot", 3, "fail to measure pixel weight");

  • trunk/psphot/src/psphotMergeSources.c

    r21251 r21519  
    122122    return detections;
    123123}
     124
     125// generate the detection structure for the supplied array of sources
     126bool psphotSetSourceParams (pmConfig *config, psArray *sources, pmPSF *psf) {
     127
     128    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     129    if (!recipe) {
     130        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
     131        return false;
     132    }
     133
     134    for (int i = 0; i < sources->n; i++) {
     135        pmSource *source = sources->data[i];
     136        pmModel *model = source->modelPSF;
     137
     138        if (source->mode & PSF_SOURCE_MASK || !isfinite(source->psfMag)) {
     139            continue;
     140        }
     141
     142        float flux = powf(10.0, -0.4 * source->psfMag);
     143        float xpos = model->params->data.F32[PM_PAR_XPOS];
     144        float ypos = model->params->data.F32[PM_PAR_YPOS];
     145
     146        pmPeak *peak = pmPeakAlloc(xpos, ypos, 1.0, PM_PEAK_LONE);
     147        peak->xf = xpos;
     148        peak->yf = ypos;
     149        peak->flux = flux; // this are being set wrong, but does it matter?
     150
     151        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
     152          peak->SN = 1.0 / source->errMag;
     153        } else {
     154          peak->SN = 0.0;
     155        }
     156
     157        source->peak = peak;
     158    }
     159
     160    psLogMsg ("psphot", 3, "%ld PSF sources loaded", sources->n);
     161
     162    return true;
     163}

  • trunk/psphot/src/psphotReplaceUnfit.c

    r21183 r21519  
    1717    replace:
    1818        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    19         source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
     19        source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    2020    }
    2121    psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%ld objects)\n", psTimerMark ("psphot.replace"), sources->n);
    … …  
    3838
    3939      // replace other sources?
    40       if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
     40      if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
    4141
    4242      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    43       source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
     43      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    4444    }
    4545    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
    … …  
    6262
    6363      // replace other sources?
    64       if (source->mode & PM_SOURCE_MODE_SUBTRACTED) continue;
     64      if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
    6565
    6666      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    67       source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     67      source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    6868    }
    6969    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
    … …  
    7575
    7676    // what is current state? (true : add; false : sub)
    77     bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
     77    bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    7878    if (state && useState) return true;
    7979
    … …  
    8686
    8787    // what is current state? (true : sub; false : add)
    88     bool state = (source->mode & PM_SOURCE_MODE_SUBTRACTED);
     88    bool state = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    8989    if (state && useState) return true;
    9090
    … …  
    9797
    9898    // what is desired state? (true : add; false : sub)
    99     bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
     99    bool newState = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    100100    if (curState == newState) return true;
    101101

  • trunk/psphot/src/psphotSourceFits.c

    r21183 r21519  
    120120        pmSourceCacheModel (blend, maskVal);
    121121        pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal);
    122         blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
     122        blend->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
     123        blend->mode |=  PM_SOURCE_MODE_BLEND_FIT;
    123124    }
    124125    NfitBlend += modelSet->n;
    … …  
    143144    pmSourceCacheModel (source, maskVal);
    144145    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    145     source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
     146    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
     147    source->mode |=  PM_SOURCE_MODE_BLEND_FIT;
    146148    return true;
    147149}
    … …  
    184186    pmSourceCacheModel (source, maskVal);
    185187    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    186     source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
     188    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
    187189    return true;
    188190}
    … …  
    285287    pmSourceCacheModel (source, maskVal);
    286288    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    287     source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
     289    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
    288290    psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
    289291    return true;
    … …  
    296298    psFree (source->modelPSF);
    297299    source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
    298     source->mode    |= PM_SOURCE_MODE_SUBTRACTED;
    299     source->mode    |= PM_SOURCE_MODE_PAIR;
     300    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     301    source->mode     |= PM_SOURCE_MODE_PAIR;
    300302
    301303    // copy most data from the primary source (modelEXT, blends stay NULL)

  • trunk/psphot/src/psphotSourceSize.c

    r21366 r21519  
    5656
    5757        // source must have been subtracted
    58         if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) {
     58        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     59            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    5960            psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n");
    6061            continue;
    … …  
    6364        psF32 **resid  = source->pixels->data.F32;
    6465        psF32 **variance = source->variance->data.F32;
    65         psImageMaskType **mask    = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     66        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    6667
    6768        // check for extendedness: measure the delta flux significance at the 1 sigma contour
    … …  
    8283        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
    8384            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
     85            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    8486            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    8587            continue;
    … …  
    98100        if (!keep) {
    99101            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
     102            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    100103            continue;
    101104        }
    … …  
    156159        }
    157160        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
    158         if (!isfinite(source->crNsigma)) continue;
     161        if (!isfinite(source->crNsigma)) {
     162            continue;
     163        }
    159164
    160165        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    … …  
    264269    // replace the source flux
    265270    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    266     source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
     271    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    267272
    268273    // flag this as a CR

  • trunk/psphot/src/psphotSourceStats.c

    r21359 r21519  
    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

  • trunk/psphot/src/psphotVisual.c

    r21366 r21519  
    11011101    Graphdata graphdata;
    11021102
    1103     bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
     1103    bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    11041104    psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    11051105
Note: See TracChangeset for help on using the changeset viewer.