IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35639


Ignore:
Timestamp:
Jun 9, 2013, 3:30:23 PM (13 years ago)
Author:
eugene
Message:

add threading for pmPSFtryFitEXT and pmPSFtryFitPSF; pass chisq conv options to pmSourceFitModel

Location:
branches/eam_branches/ipp-20130509/psModules/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_PS1_V1.c

    r35560 r35639  
    11/******************************************************************************
    2  * this file defines the PS1_V1 source shape model.  Note that these model functions are loaded
    3  * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
    4  * models use a psVector to represent the set of parameters, with the sequence used to specify
    5  * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
    6  * specifics of the model.  All models which are used as a PSF representations share a few
    7  * parameters, for which # define names are listed in pmModel.h:
     2 * this file defines the PS1_V1 source shape model.  Note that these model functions are
     3 * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
     4 * their own.  The models use a psVector to represent the set of parameters, with the
     5 * sequence used to specify the meaning of the parameter.  The meaning of the parameters
     6 * may thus vary depending on the specifics of the model.  All models which are used as a
     7 * PSF representations share a few parameters, for which # define names are listed in
     8 * pmModel.h:
    89
    910   power-law with fitted linear term

  • branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_QGAUSS.c

    r35560 r35639  
    11/******************************************************************************
    2  * this file defines the QGAUSS source shape model (XXX need a better name!).  Note that these
    3  * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
    4  * statements of their own.  The models use a psVector to represent the set of parameters, with
    5  * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
    6  * may thus vary depending on the specifics of the model.  All models which are used a PSF
    7  * representations share a few parameters, for which # define names are listed in pmModel.h:
     2 * this file defines the QGAUSS source shape model.  Note that these model functions are
     3 * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
     4 * their own.  The models use a psVector to represent the set of parameters, with the
     5 * sequence used to specify the meaning of the parameter.  The meaning of the parameters
     6 * may thus vary depending on the specifics of the model.  All models which are used as a
     7 * PSF representations share a few parameters, for which # define names are listed in
     8 * pmModel.h:
    89
    910   power-law with fitted linear term
    … …  
    1415   * PM_PAR_XPOS 2  - X center of object
    1516   * PM_PAR_YPOS 3  - Y center of object
    16    * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
    17    * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     17   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
     18   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    1819   * PM_PAR_SXY 6   - X*Y term of elliptical contour
    1920   * PM_PAR_7   7   - amplitude of the linear component (k)
    … …  
    138139# define AR_MAX 20.0
    139140# define AR_RATIO 0.99
    140 
    141141bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    142142{
    … …  
    203203          return true;
    204204      }
    205     default:
     205      default:
    206206        psAbort("invalid choice for limits");
    207207    }
    … …  
    293293    if (flux >= PAR[PM_PAR_I0]) return 1.0;
    294294
    295     // if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    296 
    297295    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    298296    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    … …  
    307305        return ( sigma * sqrt (2.0 * z) );
    308306    }
    309 
    310307    psF64 limit = flux / PAR[PM_PAR_I0];
    311308

  • branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_TRAIL.c

    r35594 r35639  
    369369    if (!isfinite(axes.theta)) return false;
    370370
    371     float size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
    372     //    if (size > psfAxes.major) { size -= psfAxes.major; }
    373     //else { size = psfAxes.major; }
     371    float size = NAN;
     372    if (!isfinite(source->moments->Mrf)) {
     373      size = axes.major;
     374    } else {
     375      size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
     376    }
    374377
    375378    float theta, peak;

  • branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtry.h

    r30044 r35639  
    9999/** fit EXT models to all possible psf sources */
    100100bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
     101bool pmPSFtryFitEXT_Threaded (psThreadJob *job);
    101102
    102103bool pmPSFtryMakePSF (bool *pGoodFit, pmPSFtry *psfTry);
    103104
    104105bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
     106bool pmPSFtryFitPSF_Threaded (psThreadJob *job);
     107
     108bool pmPSFThreads (void);
    105109
    106110/** pmPSFtryMetric()

  • branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitEXT.c

    r35560 r35639  
    4646#include "pmSourceVisual.h"
    4747
     48bool pmPSFThreads (void) {
     49
     50    psThreadTask *task = NULL;
     51
     52    task = psThreadTaskAlloc("PSF_TRY_FIT_EXT", 6);
     53    task->function = &pmPSFtryFitEXT_Threaded;
     54    psThreadTaskAdd(task);
     55    psFree(task);
     56
     57    task = psThreadTaskAlloc("PSF_TRY_FIT_PSF", 6);
     58    task->function = &pmPSFtryFitPSF_Threaded;
     59    psThreadTaskAdd(task);
     60    psFree(task);
     61
     62    return true;
     63}
     64
     65static int Next = 0;
     66
    4867// Fit an EXT model to all candidates PSF sources.
    4968// Note: this is independent of the modeled 2D variations in the PSF.
    5069bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
    51 
    52     bool status;
    5370
    5471    psTimerStart ("psf.fit");
    … …  
    6077    maskVal |= markVal;
    6178
    62     int Next = 0;
     79    Next = 0;
    6380    for (int i = 0; i < psfTry->sources->n; i++) {
    6481
    6582        pmSource *source = psfTry->sources->data[i];
    66         if (!source->moments) {
    67             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    68             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
    69             continue;
    70         }
    71         if (!source->moments->nPixels) {
    72             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
    73             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    74             continue;
    75         }
    76         // If mask object does not exist, mark the source as bad.
    77         // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
    78         // psphot to exit with a fault.
    79         if (source->maskObj == NULL) {
    80             psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
    81             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    82             continue;
    83         }
    8483
    85         source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
    86         if (source->modelEXT == NULL) {
    87             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
    88             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    89             continue;
    90         }
     84        if (!source->moments) {
     85            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     86            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
     87            continue;
     88        }
     89        if (!source->moments->nPixels) {
     90            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
     91            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     92            continue;
     93        }
     94        // If mask object does not exist, mark the source as bad.
     95        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
     96        // psphot to exit with a fault.
     97        if (source->maskObj == NULL) {
     98            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
     99            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     100            continue;
     101        }
    91102
    92         // set object mask to define valid pixels
    93         // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
    94         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
     103        source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
     104        if (source->modelEXT == NULL) {
     105            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
     106            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     107            continue;
     108        }
    95109
    96         // fit model as EXT, not PSF
    97         status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal);
     110        // do some actual work on this source
     111        psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_EXT");
     112        psArrayAdd(job->args, 1, source);
     113        psArrayAdd(job->args, 1, psfTry);
     114        psArrayAdd(job->args, 1, options);
     115       
     116        PS_ARRAY_ADD_SCALAR(job->args, i,        PS_TYPE_S32);
    98117
    99         // clear object mask to define valid pixels
    100         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
     118        PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     119        PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
    101120
    102         // exclude the poor fits
    103         if (!status) {
    104             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
    105             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    106             continue;
    107         }
    108         Next ++;
     121# if (1)
     122        if (!psThreadJobAddPending(job)) {
     123            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
     124            return false;
     125        }
     126# else
     127        if (!pmPSFtryFitEXT_Threaded(job)) {
     128            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
     129            return false;
     130        }
     131        psFree(job);
     132# endif
    109133    }
     134
     135    // wait for the threads to finish and manage results
     136    if (!psThreadPoolWait (false, true)) {
     137        psError(PS_ERR_UNKNOWN, false, "failure to model psf");
     138        return false;
     139    }
     140
     141    // we have only supplied one type of job, so we can assume the types here
     142    psThreadJob *job = NULL;
     143    while ((job = psThreadJobGetDone()) != NULL) {
     144        // we have no returned data from this operation
     145        if (job->args->n < 1) fprintf (stderr, "error with job\n");
     146        psFree(job);
     147    }
     148
    110149    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, psfTry->sources->n);
    111150    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, psfTry->sources->n);
    … …  
    118157    return true;
    119158}
     159
     160bool pmPSFtryFitEXT_Threaded (psThreadJob *job) {
     161
     162    pmSource *source =      job->args->data[0];
     163    pmPSFtry *psfTry =      job->args->data[1];
     164    pmPSFOptions *options = job->args->data[2];
     165
     166    int i = PS_SCALAR_VALUE(job->args->data[3], S32);
     167
     168    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     169    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     170
     171    // set object mask to define valid pixels
     172    // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
     173    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
     174
     175    // fit model as EXT, not PSF
     176    bool status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal);
     177
     178    // clear object mask to define valid pixels
     179    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
     180
     181    // exclude the poor fits
     182    if (!status) {
     183        psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
     184        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     185        return true;
     186    }
     187    Next ++;
     188   
     189    return true;
     190}

  • branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitPSF.c

    r34403 r35639  
    4343#include "pmSourceVisual.h"
    4444
     45static int Npsf = 0;
     46
    4547// stage 3: Refit with fixed shape parameters.  This function uses the LMM fitting, but could
    4648// be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c)
    4749bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
    4850
    49     bool status;
    50 
    5151    psTimerStart ("psf.fit");
    5252
    … …  
    5757    maskVal |= markVal;
    5858
    59     // DEBUG code: save the PSF model fit data in detail
    60 # ifdef DEBUG
    61     char filename[64];
    62     snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
    63     FILE *f = fopen (filename, "w");
    64     psAssert (f, "failed open");
    65 # endif
    66 
    67     int Npsf = 0;
     59    Npsf = 0;
    6860    for (int i = 0; i < psfTry->sources->n; i++) {
    6961
    … …  
    7769        }
    7870
    79         // set shape for this model based on PSF
     71        // set shape for this model based on PSF
    8072        psFree (source->modelPSF);
    81         source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
    82         if (source->modelPSF == NULL) {
    83             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
    84             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
    85             continue;
    86         }
    87         // PSF fit and aperture mags use different radii
    88         source->modelPSF->fitRadius = options->fitRadius;
    89         source->apRadius = options->apRadius;
    90 
    91         // set object mask to define valid pixels for PSF model fit
    92         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
    93 
    94         // fit the PSF model to the source
    95         status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
    96 
    97         // skip poor fits
    98         if (!status) {
    99             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
    100             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
    101             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
    102             continue;
    103         }
    104 
    105         // set object mask to define valid pixels for APERTURE magnitude
    106         if (options->fitRadius != options->apRadius) {
    107             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
    108             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
     73        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
     74        if (source->modelPSF == NULL) {
     75            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
     76            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
     77            return false;
    10978        }
    11079
    111         // This function calculates the psf and aperture magnitudes
    112         status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
    113         if (!status || isnan(source->apMag)) {
    114             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
    115             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
    116             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
    117             continue;
    118         }
    119 
    120         // clear object mask to define valid pixels
    121         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
    122 
    123         psfTry->fitMag->data.F32[i] = source->psfMag;
    124         psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
    125         psfTry->metricErr->data.F32[i] = source->psfMagErr;
    126 
    127         // XXX this did not work: modifies shape of psf too much
    128         // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag;
    129 
     80        // do some actual work on this source
     81        psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_PSF");
     82        psArrayAdd(job->args, 1, source);
     83        psArrayAdd(job->args, 1, psfTry);
     84        psArrayAdd(job->args, 1, options);
     85       
     86        PS_ARRAY_ADD_SCALAR(job->args, i,        PS_TYPE_S32);
     87
     88        PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     89        PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
     90
     91# if (1)
     92        if (!psThreadJobAddPending(job)) {
     93            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
     94            return false;
     95        }
     96# else
     97        if (!pmPSFtryFitPSF_Threaded(job)) {
     98            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
     99            return false;
     100        }
     101        psFree(job);
     102# endif
     103    }
     104
     105    // wait for the threads to finish and manage results
     106    if (!psThreadPoolWait (false, true)) {
     107        psError(PS_ERR_UNKNOWN, false, "failure to model psf");
     108        return false;
     109    }
     110
     111    // we have only supplied one type of job, so we can assume the types here
     112    psThreadJob *job = NULL;
     113    while ((job = psThreadJobGetDone()) != NULL) {
     114        // we have no returned data from this operation
     115        if (job->args->n < 1) fprintf (stderr, "error with job\n");
     116        psFree(job);
     117    }
     118    psfTry->psf->nPSFstars = Npsf;
     119
     120    // DEBUG code: save the PSF model fit data in detail
    130121# ifdef DEBUG
     122
     123    char filename[64];
     124    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
     125    FILE *f = fopen (filename, "w");
     126    psAssert (f, "failed open");
     127
     128    for (int i = 0; i < psfTry->sources->n; i++) {
     129
     130        // skip masked sources
     131        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
     132
     133        pmSource *source = psfTry->sources->data[i];
     134
    131135        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
    132136                 source->peak->xf, source->peak->yf,
    … …  
    136140                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY],
    137141                 source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]);
    138 # endif
    139 
    140         psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
    141         Npsf ++;
    142     }
    143     psfTry->psf->nPSFstars = Npsf;
    144 
    145 # ifdef DEBUG
     142    }
    146143    fclose (f);
    147144# endif
    … …  
    159156    return true;
    160157}
     158
     159bool pmPSFtryFitPSF_Threaded (psThreadJob *job) {
     160
     161    pmSource *source =      job->args->data[0];
     162    pmPSFtry *psfTry =      job->args->data[1];
     163    pmPSFOptions *options = job->args->data[2];
     164
     165    int i = PS_SCALAR_VALUE(job->args->data[3], S32);
     166
     167    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     168    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     169
     170    // PSF fit and aperture mags use different radii
     171    source->modelPSF->fitRadius = options->fitRadius;
     172    source->apRadius            = options->apRadius;
     173
     174    // set object mask to define valid pixels for PSF model fit
     175    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
     176
     177    // fit the PSF model to the source
     178    bool status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
     179
     180    // skip poor fits
     181    if (!status) {
     182        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
     183        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
     184        psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
     185        return true;
     186    }
     187
     188    // set object mask to define valid pixels for APERTURE magnitude
     189    if (options->fitRadius != options->apRadius) {
     190        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
     191        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
     192    }
     193
     194    // This function calculates the psf and aperture magnitudes
     195    status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
     196    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
     197
     198    if (!status || isnan(source->apMag)) {
     199        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
     200        psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
     201        return true;
     202    }
     203
     204    psfTry->fitMag->data.F32[i] = source->psfMag;
     205    psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
     206    psfTry->metricErr->data.F32[i] = source->psfMagErr;
     207
     208    psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
     209    Npsf ++;
     210    return true;
     211}

  • branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.c

    r35628 r35639  
    6262    opt->poissonErrors = true;
    6363    opt->saveCovariance = false;
     64
     65    // we default to the old algorithm
     66    opt->gainFactorMode = 0;
     67    opt->chisqConvergence = true;
    6468
    6569    return opt;
    … …  
    241245
    242246    psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
     247    myMin->gainFactorMode = options->gainFactorMode;
     248    myMin->chisqConvergence = options->chisqConvergence;
    243249
    244250    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    … …  
    274280    model->flags |= PM_MODEL_STATUS_FITTED;
    275281    if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
    276     if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
     282
     283    if (myMin->chisqConvergence) {
     284      if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
     285    } else {
     286      if (myMin->rParSigma > myMin->minTol*nParams) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
     287    }
    277288
    278289    // get the Gauss-Newton distance for fixed model parameters

  • branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.h

    r34259 r35639  
    3535    bool poissonErrors;                 ///< use poisson errors for fits?
    3636    bool saveCovariance;
     37    int gainFactorMode;
     38    bool chisqConvergence;
    3739} pmSourceFitOptions;
    3840

  • branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitSet.c

    r34403 r35639  
    568568
    569569    psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
     570    myMin->gainFactorMode = options->gainFactorMode;
     571    myMin->chisqConvergence = options->chisqConvergence;
    570572
    571573    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

  • branches/eam_branches/ipp-20130509/psModules/src/psmodules.h

    r34403 r35639  
    158158#include <pmReadoutFake.h>
    159159#include <pmPSFEnvelope.h>
     160#include <pmThreadTools.h>
    160161
    161162#endif
Note: See TracChangeset for help on using the changeset viewer.