IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21274


Ignore:
Timestamp:
Feb 3, 2009, 1:55:26 PM (17 years ago)
Author:
eugene
Message:

modify threading to call unthreaded versions of functions for nThreads <= 1 (for testing)

Location:
branches/eam_branch_20090203/psphot/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20090203/psphot/src/psphotApResid.c

    r21183 r21274  
    11# include "psphotInternal.h"
     2
     3bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal);
    24
    35# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
    … …  
    2931        nThreads = 0;
    3032    }
    31     // nThreads = 0; // XXX until testing is complete, do not thread this function
    3233
    3334    bool measureAptrend = psMetadataLookupBool (&status, recipe, "MEASURE.APTREND");
    … …  
    9495        for (int j = 0; j < cells->n; j++) {
    9596
    96             // allocate a job -- if threads are not defined, this just runs the job
    97             psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS");
    98 
    99             psArrayAdd(job->args, 1, cells->data[j]); // sources
    100             psArrayAdd(job->args, 1, psf);
    101             PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
    102             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    103 
    104             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
    105             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
    106 
    107             if (!psThreadJobAddPending(job)) {
     97            if (nThreads > 1) {
     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 {
     116                int nskip = 0;
     117                int nfail = 0;
     118
     119                if (!psphotApResidMags_Unthreaded (&nskip, &nfail, sources, psf, photMode, maskVal)) {
     120                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     121                    return false;
     122                }
     123                Nskip += nskip;
     124                Nfail += nfail;
     125            }           
     126
     127        }
     128
     129        if (nThreads > 1) {
     130            // wait for the threads to finish and manage results
     131            if (!psThreadPoolWait (false)) {
    108132                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    109                 psFree (job);
    110133                return false;
    111134            }
    112             psFree(job);
    113 
    114         }
    115 
    116         // wait for the threads to finish and manage results
    117         if (!psThreadPoolWait (false)) {
    118             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    119             return false;
    120         }
    121 
    122         // we have only supplied one type of job, so we can assume the types here
    123         psThreadJob *job = NULL;
    124         while ((job = psThreadJobGetDone()) != NULL) {
    125             if (job->args->n < 1) {
    126                 fprintf (stderr, "error with job\n");
    127             } else {
    128                 psScalar *scalar = NULL;
    129                 scalar = job->args->data[4];
    130                 Nskip += scalar->data.S32;
    131                 scalar = job->args->data[5];
    132                 Nfail += scalar->data.S32;
     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);
    133149            }
    134             psFree(job);
    135150        }
    136151    }
    … …  
    493508    return true;
    494509}
     510
     511bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
     512
     513    int Nskip = 0;
     514    int Nfail = 0;
     515
     516    for (int i = 0; i < sources->n; i++) {
     517        pmSource *source = (pmSource *) sources->data[i];
     518
     519        if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
     520        if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
     521        if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
     522        if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
     523        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
     524
     525        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
     526            Nskip ++;
     527            psTrace ("psphot", 3, "skip : bad source mag");
     528            continue;
     529        }
     530
     531        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
     532            Nfail ++;
     533            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
     534            continue;
     535        }
     536    }
     537
     538    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     539    *nskip = Nskip;
     540    *nfail = Nfail;
     541
     542    return true;
     543}

  • branches/eam_branch_20090203/psphot/src/psphotBlendFit.c

    r21250 r21274  
    11# include "psphotInternal.h"
     2
     3bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources);
    24
    35// XXX I don't like this name
    … …  
    2123        nThreads = 0;
    2224    }
    23     // nThreads = 0; // XXX until testing is complete, do not thread this function
    2425
    2526    // source fitting parameters for extended source fits
    … …  
    6364        for (int j = 0; j < cells->n; j++) {
    6465
    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)) {
     66            if (nThreads > 1) {
     67                // allocate a job -- if threads are not defined, this just runs the job
     68                psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
     69                psArray *newSources = psArrayAllocEmpty(16);
     70
     71                psArrayAdd(job->args, 1, readout);
     72                psArrayAdd(job->args, 1, recipe);
     73                psArrayAdd(job->args, 1, cells->data[j]); // sources
     74                psArrayAdd(job->args, 1, psf);
     75                psArrayAdd(job->args, 1, newSources); // return for new sources
     76                psFree (newSources);
     77
     78                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
     79                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
     80                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     81                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     82
     83                if (!psThreadJobAddPending(job)) {
     84                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     85                    psFree (job);
     86                    return NULL;
     87                }
     88                psFree(job);
     89            } else {
     90                int nfit = 0;
     91                int npsf = 0;
     92                int next = 0;
     93                int nfail = 0;
     94                psArray *newSources = psArrayAllocEmpty(16);
     95
     96                if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
     97                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     98                    return NULL;
     99                }
     100                Nfit += nfit;
     101                Npsf += npsf;
     102                Next += next;
     103                Nfail += nfail;
     104
     105                // add these back onto sources
     106                for (int k = 0; k < newSources->n; k++) {
     107                    psArrayAdd (sources, 16, newSources->data[k]);
     108                }
     109            }
     110        }
     111
     112        if (nThreads > 1) {
     113            // wait for the threads to finish and manage results
     114            if (!psThreadPoolWait (false)) {
    82115                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    83                 psFree (job);
    84116                return NULL;
    85117            }
    86             psFree(job);
    87 
    88         }
    89 
    90         // wait for the threads to finish and manage results
    91         if (!psThreadPoolWait (false)) {
    92             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    93             return NULL;
    94         }
    95 
    96         // we have only supplied one type of job, so we can assume the types here
    97         psThreadJob *job = NULL;
    98         while ((job = psThreadJobGetDone()) != NULL) {
    99             if (job->args->n < 1) {
    100                 fprintf (stderr, "error with job\n");
    101             } else {
    102                 psScalar *scalar = NULL;
    103                 scalar = job->args->data[5];
    104                 Nfit += scalar->data.S32;
    105                 scalar = job->args->data[6];
    106                 Npsf += scalar->data.S32;
    107                 scalar = job->args->data[7];
    108                 Next += scalar->data.S32;
    109                 scalar = job->args->data[8];
    110                 Nfail += scalar->data.S32;
    111 
    112                 // add these back onto sources
    113                 psArray *newSources = job->args->data[4];
    114                 for (int j = 0; j < newSources->n; j++) {
    115                     psArrayAdd (sources, 16, newSources->data[j]);
     118
     119            // we have only supplied one type of job, so we can assume the types here
     120            psThreadJob *job = NULL;
     121            while ((job = psThreadJobGetDone()) != NULL) {
     122                if (job->args->n < 1) {
     123                    fprintf (stderr, "error with job\n");
     124                } else {
     125                    psScalar *scalar = NULL;
     126                    scalar = job->args->data[5];
     127                    Nfit += scalar->data.S32;
     128                    scalar = job->args->data[6];
     129                    Npsf += scalar->data.S32;
     130                    scalar = job->args->data[7];
     131                    Next += scalar->data.S32;
     132                    scalar = job->args->data[8];
     133                    Nfail += scalar->data.S32;
     134
     135                    // add these back onto sources
     136                    psArray *newSources = job->args->data[4];
     137                    for (int j = 0; j < newSources->n; j++) {
     138                        psArrayAdd (sources, 16, newSources->data[j]);
     139                    }
    116140                }
    117             }
    118             psFree(job);
     141                psFree(job);
     142            }
    119143        }
    120144    }
    … …  
    249273    return true;
    250274}
     275
     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->mode & PM_SOURCE_MODE_SUBTRACTED) {
     339            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     340        }
     341        Nfit ++;
     342
     343        // try fitting PSFs or extended sources depending on source->mode
     344        // these functions subtract the resulting fitted source
     345        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     346            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
     347                source->type = PM_SOURCE_TYPE_EXTENDED;
     348                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
     349                Next ++;
     350                continue;
     351            }
     352        } else {
     353            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
     354                source->type = PM_SOURCE_TYPE_STAR;
     355                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
     356                Npsf ++;
     357                continue;
     358            }
     359        }
     360
     361        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
     362        Nfail ++;
     363
     364        // re-subtract the object, leave local sky
     365        pmSourceCacheModel (source, maskVal);
     366        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     367        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     368    }
     369
     370    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     371    *nfit  = Nfit;
     372    *npsf  = Npsf;
     373    *next  = Next;
     374    *nfail = Nfail;
     375
     376    return true;
     377}

  • branches/eam_branch_20090203/psphot/src/psphotGuessModels.c

    r21183 r21274  
    2222}
    2323
     24bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal);
     25
    2426// construct an initial PSF model for each object
    2527bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
    … …  
    6668        for (int j = 0; j < cells->n; j++) {
    6769
    68             // allocate a job -- if threads are not defined, this just runs the job
    69             psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
    70             psArrayAdd(job->args, 1, readout);
    71             psArrayAdd(job->args, 1, cells->data[j]); // sources
    72             psArrayAdd(job->args, 1, psf);
    73 
    74             // XXX change these to use abstract mask type info
    75             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    76             PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
    77 
    78             if (!psThreadJobAddPending(job)) {
     70            if (nThreads > 1) {
     71                // allocate a job -- if threads are not defined, this just runs the job
     72                psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
     73                psArrayAdd(job->args, 1, readout);
     74                psArrayAdd(job->args, 1, cells->data[j]); // sources
     75                psArrayAdd(job->args, 1, psf);
     76
     77                // XXX change these to use abstract mask type info
     78                PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     79                PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
     80
     81                if (!psThreadJobAddPending(job)) {
     82                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     83                    psFree (job);
     84                    return false;
     85                }
     86                psFree(job);
     87            } else {
     88                if (!psphotGuessModel_Unthreaded (readout, cells->data[j], psf, maskVal, markVal)) {
     89                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     90                    return false;
     91                }
     92            }
     93        }
     94
     95        if (nThreads > 1) {
     96            // wait for the threads to finish and manage results
     97            // wait here for the threaded jobs to finish
     98            // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
     99            if (!psThreadPoolWait (false)) {
    79100                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    80                 psFree (job);
    81101                return false;
    82102            }
    83             psFree(job);
    84         }
    85 
    86         // wait for the threads to finish and manage results
    87         // wait here for the threaded jobs to finish
    88         // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
    89         if (!psThreadPoolWait (false)) {
    90             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    91             return false;
    92         }
    93 
    94         // we have only supplied one type of job, so we can assume the types here
    95         psThreadJob *job = NULL;
    96         while ((job = psThreadJobGetDone()) != NULL) {
    97             // we have no returned data from this operation
    98             if (job->args->n < 1) {
    99                 fprintf (stderr, "error with job\n");
    100             }
    101             psFree(job);
     103
     104            // we have only supplied one type of job, so we can assume the types here
     105            psThreadJob *job = NULL;
     106            while ((job = psThreadJobGetDone()) != NULL) {
     107                // we have no returned data from this operation
     108                if (job->args->n < 1) {
     109                    fprintf (stderr, "error with job\n");
     110                }
     111                psFree(job);
     112            }
    102113        }
    103114    }
    … …  
    211222}
    212223
     224// construct models only for sources in the specified region
     225bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
     226
     227    int nSrc = 0;
     228
     229    for (int i = 0; i < sources->n; i++) {
     230        pmSource *source = sources->data[i];
     231
     232        // XXXX this is just for a test: use this to mark sources for which the model is measured
     233        // check later that all are used.
     234        source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     235
     236        // skip non-astronomical objects (very likely defects)
     237        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     238        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     239        if (!source->peak) continue;
     240
     241        nSrc ++;
     242       
     243        // XXX if a source is faint, it will not have moments measured.
     244        // it must be modelled as a PSF.  In this case, we need to use
     245        // the peak centroid to get the coordinates and get the peak flux
     246        // from the image?
     247        pmModel *modelEXT;
     248        if (!source->moments) {
     249            modelEXT = wildGuess(source, psf);
     250        } else {
     251            // use the source moments, etc to guess basic model parameters
     252            modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC
     253            if (!modelEXT) {
     254                modelEXT = wildGuess(source, psf);
     255            }
     256            // these valuse are set in pmSourceModelGuess, should this rule be in there as well?
     257            if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
     258                modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx;
     259                modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My;
     260            } else {
     261                modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf;
     262                modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf;
     263            }
     264        }
     265
     266        // set PSF parameters for this model (apply 2D shape model)
     267        pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC
     268        if (modelPSF == NULL) {
     269            psError(PSPHOT_ERR_PSF, false,
     270                    "Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
     271                    source->peak->y, source->peak->x);
     272            //
     273            // Try the centre of the image
     274            //
     275            modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
     276            modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
     277            modelPSF = pmModelFromPSF (modelEXT, psf);
     278            if (modelPSF == NULL) {
     279                psError(PSPHOT_ERR_PSF, false,
     280                        "Failed to determine PSF model at centre of image");
     281                psFree(modelEXT);
     282                return false;
     283            }
     284
     285            source->mode |= PM_SOURCE_MODE_BADPSF;
     286        }
     287        psFree (modelEXT);
     288
     289        // XXX need to define the guess flux?
     290        // set the fit radius based on the object flux limit and the model
     291        // this function affects the mask pixels
     292        psphotCheckRadiusPSF (readout, source, modelPSF, markVal);
     293
     294        // set the source PSF model
     295        source->modelPSF = modelPSF;
     296        source->modelPSF->residuals = psf->residuals;
     297
     298        pmSourceCacheModel (source, maskVal);
     299
     300    }
     301
     302    return true;
     303}

  • branches/eam_branch_20090203/psphot/src/psphotMagnitudes.c

    r21252 r21274  
    11# include "psphotInternal.h"
     2
     3bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal);
    24
    35bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
    … …  
    1719        nThreads = 0;
    1820    }
    19     // nThreads = 0; // XXX until testing is complete, do not thread this function
    2021
    2122    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
    … …  
    5859        for (int j = 0; j < cells->n; j++) {
    5960
    60             // allocate a job -- if threads are not defined, this just runs the job
    61             psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES");
    62 
    63             psArrayAdd(job->args, 1, cells->data[j]); // sources
    64             psArrayAdd(job->args, 1, psf);
    65             psArrayAdd(job->args, 1, binning);
    66             psArrayAdd(job->args, 1, backModel);
    67             psArrayAdd(job->args, 1, backStdev);
    68 
    69             PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
    70             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8); // XXX change this to use abstract mask type info
    71             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
    72 
    73             if (!psThreadJobAddPending(job)) {
     61            if (nThreads > 1) {
     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 {
     82                int nap = 0;
     83                if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) {
     84                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     85                    return false;
     86                }
     87                Nap += nap;
     88            }
     89        }
     90
     91        if (nThreads > 1) {
     92            // wait for the threads to finish and manage results
     93            if (!psThreadPoolWait (false)) {
    7494                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    75                 psFree (job);
    7695                return false;
    7796            }
    78             psFree(job);
    79 
    80         }
    81 
    82         // wait for the threads to finish and manage results
    83         if (!psThreadPoolWait (false)) {
    84             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    85             return false;
    86         }
    87 
    88         // we have only supplied one type of job, so we can assume the types here
    89         psThreadJob *job = NULL;
    90         while ((job = psThreadJobGetDone()) != NULL) {
    91             if (job->args->n < 1) {
    92                 fprintf (stderr, "error with job\n");
    93             } else {
    94                 psScalar *scalar = job->args->data[7];
    95                 Nap += scalar->data.S32;
    96             }
    97             psFree(job);
     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            }
    98109        }
    99110    }
    … …  
    116127    pmReadout *backStdev            = job->args->data[4];
    117128    pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[5],S32);
    118     psMaskType maskVal              = PS_SCALAR_VALUE(job->args->data[6],U8);
     129    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
    119130
    120131    for (int i = 0; i < sources->n; i++) {
    … …  
    150161    return true;
    151162}
     163
     164bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
     165
     166    bool status;
     167    int Nap = 0;
     168
     169    for (int i = 0; i < sources->n; i++) {
     170        pmSource *source = (pmSource *) sources->data[i];
     171        status = pmSourceMagnitudes (source, psf, photMode, maskVal);
     172        if (status && isfinite(source->apMag)) Nap ++;
     173
     174        if (backModel) {
     175            psAssert (binning, "if backModel is defined, so should binning be");
     176            source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
     177            if (isnan(source->sky) && false) {
     178                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
     179            }
     180        } else {
     181            source->sky = NAN;
     182        }
     183
     184        if (backStdev) {
     185            psAssert (binning, "if backStdev is defined, so should binning be");
     186            source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
     187            if (isnan(source->skyErr) && false) {
     188                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
     189            }
     190        } else {
     191            source->skyErr = NAN;
     192        }
     193    }
     194
     195    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     196    *nap = Nap;
     197
     198    return true;
     199}
     200
     201# if (0)
     202bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
     203
     204    bool status = false;
     205    int Nap = 0;
     206
     207    psTimerStart ("psphot.mags");
     208
     209    // select the appropriate recipe information
     210    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     211    assert (recipe);
     212
     213    // determine the number of allowed threads
     214    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     215    if (!status) {
     216        nThreads = 0;
     217    }
     218    nThreads = 0; // XXX until testing is complete, do not thread this function
     219
     220    // bit-masks to test for good/bad pixels
     221    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
     222    assert (maskVal);
     223
     224    // bit-mask to mark pixels not used in analysis
     225    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     226    assert (markVal);
     227
     228    // maskVal is used to test for rejected pixels, and must include markVal
     229    maskVal |= markVal;
     230
     231    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_APCORR | PM_SOURCE_PHOT_WEIGHT;
     232    if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH;
     233    if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP;
     234
     235    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     236    int Cx = 1, Cy = 1;
     237    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     238
     239    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     240
     241    for (int i = 0; i < cellGroups->n; i++) {
     242
     243        psArray *cells = cellGroups->data[i];
     244
     245        for (int j = 0; j < cells->n; j++) {
     246
     247            // allocate a job -- if threads are not defined, this just runs the job
     248            psThreadJob *job = psThreadJobAlloc ("PSPHOT_PSF_WEIGHTS");
     249
     250            psArrayAdd(job->args, 1, cells->data[j]); // sources
     251            psArrayAdd(job->args, 1, psf);
     252            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     253
     254            if (!psThreadJobAddPending(job)) {
     255                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     256                psFree (job);
     257                return false;
     258            }
     259            psFree(job);
     260
     261        }
     262
     263        // wait for the threads to finish and manage results
     264        if (!psThreadPoolWait (false)) {
     265            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     266            return false;
     267        }
     268
     269        // we have only supplied one type of job, so we can assume the types here
     270        psThreadJob *job = NULL;
     271        while ((job = psThreadJobGetDone()) != NULL) {
     272            if (job->args->n < 1) {
     273                fprintf (stderr, "error with job\n");
     274            }
     275            psFree(job);
     276        }
     277    }
     278
     279    psFree (cellGroups);
     280
     281    psLogMsg ("psphot", PS_LOG_DETAIL, "measure psf weights : %f sec for %ld objects\n", psTimerMark ("psphot.mags"), sources->n);
     282    return true;
     283}
     284
     285bool psphotPSFWeights_Threaded (psThreadJob *job) {
     286
     287    bool status;
     288    bool isPSF;
     289
     290    psArray *sources                = job->args->data[0];
     291    pmPSF *psf                      = job->args->data[1];
     292    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
     293
     294    for (int i = 0; i < sources->n; i++) {
     295        pmSource *source = (pmSource *) sources->data[i];
     296
     297        // we must have a valid model
     298        pmModel *model = pmSourceGetModel (&isPSF, source);
     299        if (model == NULL) {
     300          psTrace ("psphot", 3, "fail mag : no valid model");
     301          source->pixWeight = NAN;
     302          continue;
     303        }
     304
     305        status = pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal);
     306        if (!status) {
     307          psTrace ("psphot", 3, "fail to measure pixel weight");
     308          source->pixWeight = NAN;
     309          continue;
     310        }
     311
     312    }
     313
     314    return true;
     315}
     316
     317# endif
     318

  • branches/eam_branch_20090203/psphot/src/psphotSourceStats.c

    r21183 r21274  
    11# include "psphotInternal.h"
     2
     3bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe);
    24
    35psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
    … …  
    7880        for (int j = 0; j < cells->n; j++) {
    7981
    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)) {
     82            if (nThreads > 1) {
     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 {
     98                int nfail = 0;
     99                int nmoments = 0;
     100                if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) {
     101                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     102                    return NULL;
     103                }
     104                Nfail += nfail;
     105                Nmoments += nmoments;
     106            }
     107        }
     108
     109        if (nThreads > 1) {
     110            // wait for the threads to finish and manage results
     111            if (!psThreadPoolWait (false)) {
    89112                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    90                 psFree (job);
    91113                return NULL;
    92114            }
    93             psFree(job);
    94 
    95         }
    96 
    97         // wait for the threads to finish and manage results
    98         if (!psThreadPoolWait (false)) {
    99             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    100             return NULL;
    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             if (job->args->n < 1) {
    107                 fprintf (stderr, "error with job\n");
    108             } else {
    109                 psScalar *scalar = NULL;
    110                 scalar = job->args->data[2];
    111                 Nmoments += scalar->data.S32;
    112                 scalar = job->args->data[3];
    113                 Nfail += scalar->data.S32;
     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);
    114129            }
    115             psFree(job);
    116130        }
    117131    }
    … …  
    213227    return true;
    214228}
     229
     230bool 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)
     305    *nmoments = Nmoments;
     306    *nfail = Nfail;
     307   
     308    return true;
     309}
Note: See TracChangeset for help on using the changeset viewer.