IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5958


Ignore:
Timestamp:
Jan 9, 2006, 6:46:03 PM (20 years ago)
Author:
magnier
Message:

updated objects code with ApTrend concept (was in eam_rel9_b0)

Location:
branches/eam_rel9_p0/psModules/src/objects
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_rel9_p0/psModules/src/objects/Makefile.am

    r5765 r5958  
    1515        models/pmModel_QGAUSS.c \
    1616        models/pmModel_SGAUSS.c
    17    
     17
    1818psmoduleincludedir = $(includedir)
    1919psmoduleinclude_HEADERS = \

  • branches/eam_rel9_p0/psModules/src/objects/models/pmModel_PGAUSS.c

    r5257 r5958  
    2626
    2727    if (deriv != NULL) {
    28         // note difference from a pure gaussian: q = PAR[1]*r
     28        psF32 *dPAR = deriv->data.F32;
    2929        psF32 q = PAR[1]*r*r*t;
    30         deriv->data.F32[0] = +1.0;
    31         deriv->data.F32[1] = +r;
    32         deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
    33         deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
    34         deriv->data.F32[4] = -2.0*q*px*X;
    35         deriv->data.F32[5] = -2.0*q*py*Y;
    36         deriv->data.F32[6] = -q*X*Y;
     30        dPAR[0] = +1.0;
     31        dPAR[1] = +r;
     32        dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
     33        dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
     34        dPAR[4] = -2.0*q*px*X;
     35        dPAR[5] = -2.0*q*py*Y;
     36        dPAR[6] = -q*X*Y;
    3737    }
    3838    return(f);
    … …  
    4747
    4848    beta_lim[0][0].data.F32[0] = 1000;
    49     beta_lim[0][0].data.F32[1] = 10000;
     49    beta_lim[0][0].data.F32[1] = 3e6;
    5050    beta_lim[0][0].data.F32[2] = 5;
    5151    beta_lim[0][0].data.F32[3] = 5;
    … …  
    6363
    6464    params_max[0][0].data.F32[0] = 1e5;
    65     params_max[0][0].data.F32[1] = 1e6;
     65    params_max[0][0].data.F32[1] = 1e8;
    6666    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    6767    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
    … …  
    130130    params[5] = 1.2 / moments->Sy;
    131131    params[6] = 0.0;
     132
    132133    return(true);
    133134}

  • branches/eam_rel9_p0/psModules/src/objects/models/pmModel_QGAUSS.c

    r5257 r5958  
    3131
    3232    if (deriv != NULL) {
     33        psF32 *dPAR = deriv->data.F32;
     34
    3335        // note difference from a pure gaussian: q = params->data.F32[1]*r
    3436        psF32 t = PAR[1]*r*r;
    3537        psF32 q = t*(PAR[7] + 2.25*pow(z, 1.25));
    3638
    37         deriv->data.F32[0] = +1.0;
    38         deriv->data.F32[1] = +r;
    39         deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
    40         deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
    41         deriv->data.F32[4] = -2.0*q*px*X;
    42         deriv->data.F32[5] = -2.0*q*py*Y;
    43         deriv->data.F32[6] = -q*X*Y;
    44         deriv->data.F32[7] = -t*z;
     39        dPAR[0] = +1.0;
     40        dPAR[1] = +r;
     41        dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
     42        dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
     43        dPAR[4] = -2.0*q*px*X;
     44        dPAR[5] = -2.0*q*py*Y;
     45        dPAR[6] = -q*X*Y;
     46        dPAR[7] = -t*z;
    4547    }
    4648    return(f);
    … …  
    5557
    5658    beta_lim[0][0].data.F32[0] = 1000;
    57     beta_lim[0][0].data.F32[1] = 10000;
     59    beta_lim[0][0].data.F32[1] = 3e6;
    5860    beta_lim[0][0].data.F32[2] = 5;
    5961    beta_lim[0][0].data.F32[3] = 5;
    … …  
    7375
    7476    params_max[0][0].data.F32[0] = 1e5;
    75     params_max[0][0].data.F32[1] = 1e6;
     77    params_max[0][0].data.F32[1] = 1e8;
    7678    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    7779    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
    … …  
    198200    status &= ((dPAR[1]/PAR[1]) < 0.5);
    199201
    200     if (status)
    201         return true;
    202     return false;
    203 }
     202    if (!status)
     203        return false;
     204    return true;
     205}

  • branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.c

    r5255 r5958  
    1010#include "models/pmModel_SGAUSS.c"
    1111
    12 static pmModelGroup models[] = {
    13                                    {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
    14                                    {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
    15                                    {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
    16                                    {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
    17                                };
     12static pmModelGroup defaultModels[] = {
     13                                          {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
     14                                          {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
     15                                          {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
     16                                          {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
     17                                      };
     18
     19static pmModelGroup *models = NULL;
     20static int Nmodels = 0;
     21
     22static void ModelGroupFree (pmModelGroup *modelGroup)
     23{
     24
     25    if (modelGroup == NULL)
     26        return;
     27    psFree (modelGroup);
     28    return;
     29}
     30
     31pmModelGroup *pmModelGroupAlloc (int nModels)
     32{
     33
     34    pmModelGroup *modelGroup = (pmModelGroup *) psAlloc (nModels * sizeof(pmModelGroup));
     35    psMemSetDeallocator(modelGroup, (psFreeFunc) ModelGroupFree);
     36    return (modelGroup);
     37}
     38
     39void pmModelGroupAdd (pmModelGroup *model)
     40{
     41
     42    if (models == NULL) {
     43        pmModelGroupInit ();
     44    }
     45
     46    Nmodels ++;
     47    models = (pmModelGroup *) psRealloc (models, Nmodels*sizeof(pmModelGroup));
     48    models[Nmodels-1] = model[0];
     49    return;
     50}
     51
     52void pmModelGroupInit (void)
     53{
     54
     55    int Nnew = sizeof (defaultModels) / sizeof (pmModelGroup);
     56
     57    models = pmModelGroupAlloc (Nnew);
     58    for (int i = 0; i < Nnew; i++) {
     59        models[i] = defaultModels[i];
     60    }
     61    Nmodels = Nnew;
     62    return;
     63}
    1864
    1965pmModelFunc pmModelFunc_GetFunction (pmModelType type)
    2066{
    21     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    2267    if ((type < 0) || (type >= Nmodels)) {
    2368        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    2974pmModelFlux pmModelFlux_GetFunction (pmModelType type)
    3075{
    31     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    3276    if ((type < 0) || (type >= Nmodels)) {
    3377        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    3983pmModelRadius pmModelRadius_GetFunction (pmModelType type)
    4084{
    41     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    4285    if ((type < 0) || (type >= Nmodels)) {
    4386        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    4992pmModelLimits pmModelLimits_GetFunction (pmModelType type)
    5093{
    51     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    5294    if ((type < 0) || (type >= Nmodels)) {
    5395        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    59101pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type)
    60102{
    61     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    62103    if ((type < 0) || (type >= Nmodels)) {
    63104        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    69110pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type)
    70111{
    71     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    72112    if ((type < 0) || (type >= Nmodels)) {
    73113        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    79119pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type)
    80120{
    81     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    82121    if ((type < 0) || (type >= Nmodels)) {
    83122        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    89128psS32 pmModelParameterCount (pmModelType type)
    90129{
    91     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    92130    if ((type < 0) || (type >= Nmodels)) {
    93131        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    … …  
    99137psS32 pmModelSetType (char *name)
    100138{
    101 
    102     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    103139    for (int i = 0; i < Nmodels; i++) {
    104140        if (!strcmp(models[i].name, name)) {
    … …  
    111147char *pmModelGetType (pmModelType type)
    112148{
    113 
    114     int Nmodels = sizeof (models) / sizeof (pmModelGroup);
    115149    if ((type < 0) || (type >= Nmodels)) {
    116150        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");

  • branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.h

    r5255 r5958  
    99 *  @author EAM, IfA
    1010 *
    11  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-10-10 19:53:40 $
     11 *  @version $Revision: 1.1.16.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2006-01-10 04:46:03 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    209209pmModelGroup;
    210210
     211// allocate a pmModelGroup to hold nModels entries
     212pmModelGroup *pmModelGroupAlloc (int nModels);
     213
     214// initialize the internal (static) model group with the default models
     215void pmModelGroupInit (void);
     216
     217// add a new model to the internal (static) model group
     218void pmModelGroupAdd (pmModelGroup *model);
     219
    211220# endif

  • branches/eam_rel9_p0/psModules/src/objects/pmObjects.c

    r5765 r5958  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-12 21:14:38 $
     8 *  @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-01-10 04:46:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    271271    tmp->chisq = 0.0;
    272272    tmp->nIter = 0;
     273    tmp->radius = 0;
     274    tmp->status = PM_MODEL_UNTRIED;
     275
    273276    psS32 Nparams = pmModelParameterCount(type);
    274277    if (Nparams == 0) {
    … …  
    307310    tmp->mask = NULL;
    308311    tmp->moments = NULL;
     312    tmp->blends = NULL;
    309313    tmp->modelPSF = NULL;
    310314    tmp->modelFLT = NULL;
    311     tmp->type = 0;
     315    tmp->type = PM_SOURCE_UNKNOWN;
     316    tmp->mode = PM_SOURCE_DEFAULT;
    312317    psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
    313318
    … …  
    827832    psS32 numPixels = 0;
    828833    psF32 Sum = 0.0;
     834    psF32 Var = 0.0;
    829835    psF32 X1 = 0.0;
    830836    psF32 Y1 = 0.0;
    … …  
    850856    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    851857        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    852             if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
     858            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
    853859                continue;
     860            }
    854861
    855862            psF32 xDiff = col + source->pixels->col0 - xPeak;
    … …  
    858865            // XXX EAM : calculate xDiff, yDiff up front;
    859866            //           radius is just a function of (xDiff, yDiff)
    860             if (!VALID_RADIUS(xDiff, yDiff, R2))
     867            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
    861868                continue;
     869            }
    862870
    863871            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     872            psF32 wDiff = source->weight->data.F32[row][col];
    864873
    865874            // XXX EAM : check for valid S/N in pixel
    866875            // XXX EAM : should this limit be user-defined?
    867             if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1)
     876            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
    868877                continue;
    869 
     878            }
     879
     880            Var += wDiff;
    870881            Sum += pDiff;
    871882            X1  += xDiff * pDiff;
    … …  
    880891        }
    881892    }
     893
     894    // if we have less than (1/4) of the possible pixels, force a retry
    882895    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    883     if ((numPixels < 3) || (Sum <= 0)) {
    884         psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
     896    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
     897        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
    885898        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    886899        return (false);
    … …  
    899912    y = Y1/Sum;
    900913    if ((fabs(x) > radius) || (fabs(y) > radius)) {
    901         psTrace (".psModules.pmSourceMoments", 5,
     914        psTrace (".psModules.pmSourceMoments", 3,
    902915                 "large centroid swing; invalid peak %d, %d\n",
    903916                 source->peak->x, source->peak->y);
    … …  
    912925    source->moments->Sxy = XY/Sum - x*y;
    913926    source->moments->Sum = Sum;
     927    source->moments->SN  = Sum / sqrt(Var);
    914928    source->moments->Peak = peakPixel;
    915929    source->moments->nPixels = numPixels;
    … …  
    9941008        psImage *splane = NULL;
    9951009        int binX, binY;
     1010        bool status;
     1011
     1012        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
     1013        if (!status)
     1014            SX_MAX = 10.0;
     1015        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
     1016        if (!status)
     1017            SY_MAX = 10.0;
    9961018
    9971019        // construct a sigma-plane image
    9981020        // psImageAlloc does zero the data
    999         splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
     1021        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
    10001022        for (int i = 0; i < splane->numRows; i++)
    10011023        {
    … …  
    11311153XXX: How can this function ever return FALSE?
    11321154 
    1133 XXX EAM : add the saturated mask value to metadata
     1155EAM: I moved S/N calculation to pmSourceMoments, using weight image
    11341156*****************************************************************************/
    11351157
    … …  
    11461168    int Ncr      = 0;
    11471169    int Nsatstar = 0;
    1148 
     1170    // psRegion allArray = psRegionSet (0, 0, 0, 0);
     1171    psRegion inner;
     1172
     1173    // report stats on S/N values for star-like objects
    11491174    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    11501175    starsn->n = 0;
    … …  
    11521177    // check return status value (do these exist?)
    11531178    bool status;
    1154     psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
    1155     psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
    11561179    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
    1157     // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
    11581180
    11591181    // XXX allow clump size to be scaled relative to sigmas?
    … …  
    11681190        psF32 sigY = tmpSrc->moments->Sy;
    11691191
    1170         // calculate and save signal-to-noise estimates
    1171         psF32 S  = tmpSrc->moments->Sum;
    1172         psF32 A  = 4 * M_PI * sigX * sigY;
    1173         psF32 B  = tmpSrc->moments->Sky;
    1174         psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
    1175         psF32 SN = (S * sqrt(GAIN) / RT);
    1176         tmpSrc->moments->SN = SN;
    1177 
    11781192        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1179         //
    1180         // XXX: Must verify this region (the region argument was added to psImageCountPixelMask()
    1181         // after EAM wrote this code.
    1182         //
    1183         psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows);
    1184         int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED);
     1193        inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
     1194        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
    11851195
    11861196        // saturated star (size consistent with PSF or larger)
    11871197        // Nsigma should be user-configured parameter
    11881198        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
     1199        big = true;
    11891200        if ((Nsatpix > 1) && big) {
    1190             tmpSrc->type = PM_SOURCE_SATSTAR;
     1201            tmpSrc->type = PM_SOURCE_STAR;
     1202            tmpSrc->mode = PM_SOURCE_SATSTAR;
    11911203            Nsatstar ++;
    11921204            continue;
    … …  
    11961208        if (Nsatpix > 1) {
    11971209            tmpSrc->type = PM_SOURCE_SATURATED;
     1210            tmpSrc->mode = PM_SOURCE_DEFAULT;
    11981211            Nsat ++;
    11991212            continue;
    … …  
    12051218        if ((sigX < 0.05) || (sigY < 0.05)) {
    12061219            tmpSrc->type = PM_SOURCE_DEFECT;
     1220            tmpSrc->mode = PM_SOURCE_DEFAULT;
    12071221            Ncr ++;
    12081222            continue;
    … …  
    12121226        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    12131227            tmpSrc->type = PM_SOURCE_GALAXY;
     1228            tmpSrc->mode = PM_SOURCE_DEFAULT;
    12141229            Ngal ++;
    12151230            continue;
    … …  
    12171232
    12181233        // the rest are probable stellar objects
    1219         starsn->data.F32[starsn->n] = SN;
     1234        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
    12201235        starsn->n ++;
    12211236        Nstar ++;
    12221237
    1223         // PSF star (within 1.5 sigma of clump center
     1238        // PSF star (within 1.5 sigma of clump center, S/N > limit)
    12241239        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    1225         if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
    1226             tmpSrc->type = PM_SOURCE_PSFSTAR;
     1240        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
     1241            tmpSrc->type = PM_SOURCE_STAR;
     1242            tmpSrc->mode = PM_SOURCE_PSFSTAR;
    12271243            Npsf ++;
    12281244            continue;
    … …  
    12301246
    12311247        // random type of star
    1232         tmpSrc->type = PM_SOURCE_OTHER;
     1248        tmpSrc->type = PM_SOURCE_STAR;
     1249        tmpSrc->mode = PM_SOURCE_DEFAULT;
    12331250    }
    12341251
    … …  
    15511568#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
    15521569#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    1553 /******************************************************************************
    1554 pmSourceFitModel(source, model): must create the appropiate arguments to the
    1555 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
    1556  
    1557 XXX: should there be a mask value?
    1558 XXX EAM : fit the specified model (not necessarily the one in source)
    1559 *****************************************************************************/
    1560 bool pmSourceFitModel_v5(pmSource *source,
    1561                          pmModel *model,
    1562                          const bool PSF)
     1570bool pmSourceFitModel (pmSource *source,
     1571                       pmModel *model,
     1572                       const bool PSF)
    15631573{
    15641574    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    … …  
    15691579    PS_ASSERT_PTR_NON_NULL(source->mask, false);
    15701580    PS_ASSERT_PTR_NON_NULL(source->weight, false);
     1581
     1582    // XXX EAM : is it necessary for the mask & weight to exist?  the
     1583    //           tests below could be conditions (!NULL)
     1584
    15711585    psBool fitStatus = true;
    15721586    psBool onPic     = true;
    15731587    psBool rc        = true;
    15741588
    1575     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1576     //           tests below could be conditions (!NULL)
    1577 
    15781589    psVector *params = model->params;
    15791590    psVector *dparams = model->dparams;
    … …  
    15821593    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    15831594
    1584     int nParams = PSF ? params->n - 4 : params->n;
    1585 
    1586     // find the number of valid pixels
    1587     // XXX EAM : this loop and the loop below could just be one pass
    1588     //           using the psArrayAdd and psVectorExtend functions
    1589     psS32 count = 0;
     1595    int nParams = PSF ? 4 : params->n;
     1596
     1597    // maximum number of valid pixels
     1598    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
     1599
     1600    // construct the coordinate and value entries
     1601    psArray *x = psArrayAlloc(nPix);
     1602    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
     1603    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
     1604
     1605    nPix = 0;
    15901606    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    15911607        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1592             if (source->mask->data.U8[i][j] == 0) {
    1593                 count++;
    1594             }
    1595         }
    1596     }
    1597     if (count <  nParams + 1) {
     1608            if (source->mask->data.U8[i][j]) {
     1609                continue;
     1610            }
     1611            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1612
     1613            // Convert i/j to image space:
     1614            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1615            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     1616            x->data[nPix] = (psPtr *) coord;
     1617            y->data.F32[nPix] = source->pixels->data.F32[i][j];
     1618
     1619            // psMinimizeLMChi2 takes wt = 1/dY^2
     1620            if (source->weight->data.F32[i][j] == 0) {
     1621                continue;
     1622            }
     1623            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
     1624            nPix++;
     1625        }
     1626    }
     1627    if (nPix <  nParams + 1) {
    15981628        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    15991629        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     1630        model->status = PM_MODEL_BADARGS;
     1631        psFree (x);
     1632        psFree (y);
     1633        psFree (yErr);
    16001634        return(false);
    16011635    }
    1602 
    1603     // construct the coordinate and value entries
    1604     psArray *x = psArrayAlloc(count);
    1605     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    1606     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    1607     psS32 tmpCnt = 0;
    1608     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1609         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1610             if (source->mask->data.U8[i][j] == 0) {
    1611                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1612                 // XXX: Convert i/j to image space:
    1613                 // XXX EAM: coord order is (x,y) == (col,row)
    1614                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1615                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1616                 x->data[tmpCnt] = (psPtr *) coord;
    1617                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1618                 yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
    1619                 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
    1620                 //           the minimization function calculates sq()
    1621                 tmpCnt++;
    1622             }
    1623         }
    1624     }
     1636    x->n = nPix;
     1637    y->n = nPix;
     1638    yErr->n = nPix;
    16251639
    16261640    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    … …  
    16381652    }
    16391653
    1640     // XXX EAM : covar must be F64?
    1641     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    1642 
    1643     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    1644     fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
    1645     for (int i = 0; i < dparams->n; i++) {
    1646         if ((paramMask != NULL) && paramMask->data.U8[i])
    1647             continue;
    1648         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    1649     }
    1650 
    1651     // XXX EAM: we need to do something (give an error?) if rc is false
    1652     // XXX EAM: psMinimizeLMChi2 does not check convergence
    1653 
    1654     // XXX models can go insane: reject these
    1655     onPic &= (params->data.F32[2] >= source->pixels->col0);
    1656     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1657     onPic &= (params->data.F32[3] >= source->pixels->row0);
    1658     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1659 
    1660     // XXX EAM: save the resulting chisq, nDOF, nIter
    1661     model->chisq = myMin->value;
    1662     model->nIter = myMin->iter;
    1663     model->nDOF  = y->n - nParams;
    1664 
    1665     // XXX EAM get the Gauss-Newton distance for fixed model parameters
    1666     if (paramMask != NULL) {
    1667         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    1668         psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
    1669         for (int i = 0; i < dparams->n; i++) {
    1670             if (!paramMask->data.U8[i])
    1671                 continue;
    1672             dparams->data.F32[i] = delta->data.F64[i];
    1673         }
    1674         psFree (delta);
    1675     }
    1676 
    1677     psFree(x);
    1678     psFree(y);
    1679     psFree(yErr);
    1680     psFree(myMin);
    1681     psFree(covar);
    1682     psFree(paramMask);
    1683 
    1684     rc = (onPic && fitStatus);
    1685     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1686     return(rc);
    1687 }
    1688 
    1689 // XXX EAM : new version with parameter range limits and weight enhancement
    1690 bool pmSourceFitModel (pmSource *source,
    1691                        pmModel *model,
    1692                        const bool PSF)
    1693 {
    1694     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1695     PS_ASSERT_PTR_NON_NULL(source, false);
    1696     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1697     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1698     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1699     PS_ASSERT_PTR_NON_NULL(source->mask, false);
    1700     PS_ASSERT_PTR_NON_NULL(source->weight, false);
    1701 
    1702     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1703     //           tests below could be conditions (!NULL)
    1704 
    1705     psBool fitStatus = true;
    1706     psBool onPic     = true;
    1707     psBool rc        = true;
    1708     psF32  Ro, ymodel;
    1709 
    1710     psVector *params = model->params;
    1711     psVector *dparams = model->dparams;
    1712     psVector *paramMask = NULL;
    1713 
    1714     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    1715 
    1716     // XXX EAM : I need to use the sky value to constrain the weight model
    1717     int nParams = PSF ? params->n - 4 : params->n;
    1718     psF32 So = params->data.F32[0];
    1719 
    1720     // find the number of valid pixels
    1721     // XXX EAM : this loop and the loop below could just be one pass
    1722     //           using the psArrayAdd and psVectorExtend functions
    1723     psS32 count = 0;
    1724     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1725         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1726             if (source->mask->data.U8[i][j] == 0) {
    1727                 count++;
    1728             }
    1729         }
    1730     }
    1731     if (count <  nParams + 1) {
    1732         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    1733         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    1734         return(false);
    1735     }
    1736 
    1737     // construct the coordinate and value entries
    1738     psArray *x = psArrayAlloc(count);
    1739     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    1740     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    1741     psS32 tmpCnt = 0;
    1742     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1743         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1744             if (source->mask->data.U8[i][j] == 0) {
    1745                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1746                 // XXX: Convert i/j to image space:
    1747                 // XXX EAM: coord order is (x,y) == (col,row)
    1748                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1749                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1750                 x->data[tmpCnt] = (psPtr *) coord;
    1751                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1752 
    1753                 // compare observed flux to model flux to adjust weight
    1754                 ymodel = modelFunc (NULL, model->params, coord);
    1755 
    1756                 // this test enhances the weight based on deviation from the model flux
    1757                 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
    1758 
    1759                 // psMinimizeLMChi2_EAM takes wt = 1/dY^2
    1760                 if (source->weight->data.F32[i][j] == 0) {
    1761                     yErr->data.F32[tmpCnt] = 0.0;
    1762                 } else {
    1763                     yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
    1764                 }
    1765                 tmpCnt++;
    1766             }
    1767         }
    1768     }
    1769 
    1770     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    1771                             PM_SOURCE_FIT_MODEL_TOLERANCE);
    1772 
    1773     // PSF model only fits first 4 parameters, FLT model fits all
    1774     if (PSF) {
    1775         paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    1776         for (int i = 0; i < 4; i++) {
    1777             paramMask->data.U8[i] = 0;
    1778         }
    1779         for (int i = 4; i < paramMask->n; i++) {
    1780             paramMask->data.U8[i] = 1;
    1781         }
    1782     }
    1783 
    1784     // XXX EAM : I've added three types of parameter range checks
    1785     // XXX EAM : this requires my new psMinimization functions
     1654    // Set the parameter range checks
    17861655    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    17871656    psVector *beta_lim = NULL;
    … …  
    18071676    }
    18081677
    1809     // XXX EAM: we need to do something (give an error?) if rc is false
    1810     // XXX EAM: psMinimizeLMChi2 does not check convergence
    1811 
    1812     // XXX models can go insane: reject these
    1813     onPic &= (params->data.F32[2] >= source->pixels->col0);
    1814     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1815     onPic &= (params->data.F32[3] >= source->pixels->row0);
    1816     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1817 
    1818     // XXX EAM: save the resulting chisq, nDOF, nIter
     1678    // save the resulting chisq, nDOF, nIter
    18191679    model->chisq = myMin->value;
    18201680    model->nIter = myMin->iter;
    18211681    model->nDOF  = y->n - nParams;
    18221682
    1823     // XXX EAM get the Gauss-Newton distance for fixed model parameters
     1683    // get the Gauss-Newton distance for fixed model parameters
    18241684    if (paramMask != NULL) {
    18251685        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    … …  
    18321692    }
    18331693
    1834     psFree(paramMask);
     1694    // set the model success or failure status
     1695    if (!fitStatus) {
     1696        model->status = PM_MODEL_NONCONVERGE;
     1697    } else {
     1698        model->status = PM_MODEL_SUCCESS;
     1699    }
     1700
     1701    // models can go insane: reject these
     1702    onPic &= (params->data.F32[2] >= source->pixels->col0);
     1703    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
     1704    onPic &= (params->data.F32[3] >= source->pixels->row0);
     1705    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     1706    if (!onPic) {
     1707        model->status = PM_MODEL_OFFIMAGE;
     1708    }
     1709
     1710    source->mode |= PM_SOURCE_FITTED;
     1711
    18351712    psFree(x);
    18361713    psFree(y);
     1714    psFree(yErr);
    18371715    psFree(myMin);
     1716    psFree(covar);
     1717    psFree(paramMask);
    18381718
    18391719    rc = (onPic && fitStatus);
    … …  
    18461726                             pmModel *model,
    18471727                             bool center,
    1848                              psS32 flag)
     1728                             bool sky,
     1729                             bool add
     1730                                )
    18491731{
    18501732    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    … …  
    18591741    psS32 imageCol;
    18601742    psS32 imageRow;
     1743    psF32 skyValue = params->data.F32[0];
     1744    psF32 pixelValue;
    18611745
    18621746    for (psS32 i = 0; i < image->numRows; i++) {
    … …  
    18641748            if ((mask != NULL) && mask->data.U8[i][j])
    18651749                continue;
    1866             psF32 pixelValue;
     1750
    18671751            // XXX: Should you be adding the pixels for the entire subImage,
    18681752            // or a radius of pixels around it?
    … …  
    18821766            x->data.F32[0] = (float) imageCol;
    18831767            x->data.F32[1] = (float) imageRow;
    1884             pixelValue = modelFunc (NULL, params, x);
    1885             // fprintf (stderr, "%f %f  %d %d  %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue);
    1886 
    1887             if (flag == 1) {
    1888                 pixelValue = -pixelValue;
    1889             }
    1890 
    1891             // XXX: Must figure out how to calculate the image coordinates and
    1892             // how to use the boolean "center" flag.
    1893 
    1894             image->data.F32[i][j]+= pixelValue;
     1768
     1769            // set the appropriate pixel value for this coordinate
     1770            if (sky) {
     1771                pixelValue = modelFunc (NULL, params, x);
     1772            } else {
     1773                pixelValue = modelFunc (NULL, params, x) - skyValue;
     1774            }
     1775
     1776
     1777            // add or subtract the value
     1778            if (add
     1779               ) {
     1780                image->data.F32[i][j] += pixelValue;
     1781            }
     1782            else {
     1783                image->data.F32[i][j] -= pixelValue;
     1784            }
    18951785        }
    18961786    }
    … …  
    19071797                      psImage *mask,
    19081798                      pmModel *model,
    1909                       bool center)
    1910 {
    1911     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1912     psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
     1799                      bool center,
     1800                      bool sky)
     1801{
     1802    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1803    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
    19131804    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    19141805    return(rc);
    … …  
    19201811                      psImage *mask,
    19211812                      pmModel *model,
    1922                       bool center)
    1923 {
    1924     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1925     psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
     1813                      bool center,
     1814                      bool sky)
     1815{
     1816    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1817    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
    19261818    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    19271819    return(rc);

  • branches/eam_rel9_p0/psModules/src/objects/pmObjects.h

    r5765 r5958  
    1010 *  @author GLG, MHPCC
    1111 *
    12  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-12-12 21:14:38 $
     12 *  @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-01-10 04:46:03 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    118118pmPSFClump;
    119119
     120// type of model carried by the pmModel structure
    120121typedef int pmModelType;
    121 #define PS_MODEL_GAUSS 0
    122 #define PS_MODEL_PGAUSS 1
    123 #define PS_MODEL_QGAUSS 2
    124 #define PS_MODEL_SGAUSS 3
    125 
     122
     123typedef enum {
     124    PM_MODEL_UNTRIED,               ///< model fit not yet attempted
     125    PM_MODEL_SUCCESS,               ///< model fit succeeded
     126    PM_MODEL_NONCONVERGE,           ///< model fit did not converge
     127    PM_MODEL_OFFIMAGE,              ///< model fit drove out of range
     128    PM_MODEL_BADARGS                ///< model fit called with invalid args
     129} pmModelStatus;
    126130
    127131/** pmModel data structure
    … …  
    142146    int nDOF;    ///< number of degrees of freedom
    143147    int nIter;    ///< number of iterations to reach min
     148    int status;         ///< fit status
    144149    float radius;   ///< fit radius actually used
    145150}
    … …  
    157162 */
    158163typedef enum {
     164    PM_SOURCE_UNKNOWN,                  ///< a cosmic-ray
    159165    PM_SOURCE_DEFECT,                   ///< a cosmic-ray
    160166    PM_SOURCE_SATURATED,                ///< random saturated pixels
    161 
    162     PM_SOURCE_SATSTAR,                  ///< a saturated star
    163     PM_SOURCE_PSFSTAR,                  ///< a PSF star
    164     PM_SOURCE_GOODSTAR,                 ///< a good-quality star
    165 
    166     PM_SOURCE_POOR_FIT_PSF,             ///< poor quality PSF fit
    167     PM_SOURCE_FAIL_FIT_PSF,             ///< failed to get a good PSF fit
    168     PM_SOURCE_FAINTSTAR,                ///< below S/N cutoff
    169 
     167    PM_SOURCE_STAR,                     ///< a good-quality star
    170168    PM_SOURCE_GALAXY,                   ///< an extended object (galaxy)
    171     PM_SOURCE_FAINT_GALAXY,             ///< a galaxy below S/N cutoff
    172     PM_SOURCE_DROP_GALAXY,              ///< ?
    173     PM_SOURCE_FAIL_FIT_GAL,             ///< failed on the galaxy fit
    174     PM_SOURCE_POOR_FIT_GAL,             ///< poor quality galaxy fit
    175 
    176     PM_SOURCE_OTHER,                    ///< unidentified
    177169} pmSourceType;
     170
     171typedef enum {
     172    PM_SOURCE_DEFAULT    = 0x0000, ///<
     173    PM_SOURCE_PSFMODEL   = 0x0001, ///<
     174    PM_SOURCE_FLTMODEL   = 0x0002, ///<
     175    PM_SOURCE_SUBTRACTED = 0x0004, ///<
     176    PM_SOURCE_FITTED     = 0x0008, ///<
     177    PM_SOURCE_FAIL       = 0x0010, ///<
     178    PM_SOURCE_POOR       = 0x0020, ///<
     179    PM_SOURCE_PAIR       = 0x0040, ///<
     180    PM_SOURCE_PSFSTAR    = 0x0080, ///<
     181    PM_SOURCE_SATSTAR    = 0x0100, ///<
     182    PM_SOURCE_BLEND      = 0x0200, ///<
     183    PM_SOURCE_LINEAR     = 0x0400, ///<
     184    PM_SOURCE_TEMPSUB    = 0x0800, ///< XXX get me a better name!
     185} pmSourceMode;
    178186
    179187/** pmSource data structure
    … …  
    194202    pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
    195203    pmSourceType type;   ///< Best identification of object.
     204    pmSourceMode mode;   ///< Best identification of object.
     205    psArray *blends;
     206    float apMag;
     207    float fitMag;
    196208}
    197209pmSource;
    … …  
    507519    psImage *mask,   ///< The image pixel mask (valid == 0)
    508520    pmModel *model,   ///< The input pmModel
    509     bool center    ///< A boolean flag that determines whether pixels are centered
     521    bool center,    ///< A boolean flag that determines whether pixels are centered
     522    bool sky        ///< A boolean flag that determines if the sky is subtracted
    510523);
    511524
    … …  
    525538    psImage *mask,   ///< The image pixel mask (valid == 0)
    526539    pmModel *model,   ///< The input pmModel
    527     bool center    ///< A boolean flag that determines whether pixels are centered
     540    bool center,    ///< A boolean flag that determines whether pixels are centered
     541    bool sky        ///< A boolean flag that determines if the sky is subtracted
    528542);
    529543

  • branches/eam_rel9_p0/psModules/src/objects/pmPSF.c

    r5765 r5958  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-12 21:14:38 $
     8 *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-01-10 04:46:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6969 X-center
    7070 Y-center
    71  ???: Sky background value?
    72  ???: Zo?
     71 Sky background value
     72 Object Normalization
    7373 *****************************************************************************/
    7474pmPSF *pmPSFAlloc (pmModelType type)
    … …  
    8383    psf->dApResid = 0.0;
    8484    psf->skyBias  = 0.0;
     85
     86    // the ApTrend components are (x, y, rflux)
     87    psf->ApTrend = psPolynomial3DAlloc (2, 2, 1, PS_POLYNOMIAL_ORD);
     88
     89    // we do not allow any rflux vs X,Y cross-terms
     90    psf->ApTrend->mask[0][1][1] = 1;  // rflux x^0 y^1
     91    psf->ApTrend->mask[0][2][1] = 1;  // rflux x^0 y^2
     92    psf->ApTrend->mask[1][0][1] = 1;  // rflux x^1 y^0
     93    psf->ApTrend->mask[1][1][1] = 1;  // rflux x^1 y^1
     94    psf->ApTrend->mask[1][2][1] = 1;  // rflux x^1 y^2
     95    psf->ApTrend->mask[2][0][1] = 1;  // rflux x^2 y^0
     96    psf->ApTrend->mask[2][1][1] = 1;  // rflux x^2 y^1
     97    psf->ApTrend->mask[2][2][1] = 1;  // rflux x^2 y^2
     98
     99    // only 2nd order terms, no such combinations
     100    psf->ApTrend->mask[2][2][0] = 1;  // x^2 y^2
     101    psf->ApTrend->mask[2][1][0] = 1;  // x^2 y^1
     102    psf->ApTrend->mask[1][2][0] = 1;  // x^1 y^2
     103
     104    // total free parameters = 18 - 11 = 7
    85105
    86106    Nparams = pmModelParameterCount (type);

  • branches/eam_rel9_p0/psModules/src/objects/pmPSF.h

    r5255 r5958  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-10-10 19:53:40 $
     8 *  @version $Revision: 1.1.18.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-01-10 04:46:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3333    pmModelType type;   ///< PSF Model in use
    3434    psArray *params;   ///< Model parameters (psPolynomial2D)
     35    psPolynomial3D *ApTrend;  ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst)
     36    float ApResid;   ///< ???
     37    float dApResid;   ///< ???
     38    float skyBias;   ///< ???
    3539    float chisq;   ///< PSF goodness statistic
    36     float ApResid;                      ///< ???
    37     float dApResid;                     ///< ???
    38     float skyBias;                      ///< ???
    3940    int nPSFstars;   ///< number of stars used to measure PSF
     41    int nApResid;   ///< number of stars used to measure ApResid
    4042}
    4143pmPSF;

  • branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.c

    r5765 r5958  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-12-12 21:14:38 $
     7 *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-01-10 04:46:03 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    119119    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    120120
    121     // make this optional?
    122     // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
    123 
    124121    // stage 2: construct a psf (pmPSF) from this collection of model fits
    125122    pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
    … …  
    169166
    170167    }
     168    psfTry->psf->nPSFstars = Npsf;
     169
    171170    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    172171    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
    … …  
    177176    // XXX this function wants aperture radius for pmSourcePhotometry
    178177    pmPSFtryMetric (psfTry, RADIUS);
    179     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    180               modelName,
    181               psfTry->psf->ApResid,
    182               psfTry->psf->dApResid,
    183               psfTry->psf->skyBias);
    184 
     178    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
     179              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    185180    return (psfTry);
    186181}
    187182
    188 
    189183bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
    190184{
    191 
    192     float dBin;
    193     int   nKeep, nSkip;
    194 
    195185    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    196186    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    … …  
    210200    }
    211201
    212     // find min and max of (1/flux):
    213     psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    214     psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
    215 
    216     // build binned versions of rflux, metric
    217     dBin = (stats->max - stats->min) / 10.0;
    218     psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
    219     psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
    220     psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
    221     psFree (stats);
    222 
    223     psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
    224 
    225     // group data in daBin bins, measure lower 50% mean
    226     for (int i = 0; i < daBin->n; i++) {
    227 
    228         psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    229         tmp->n = 0;
    230 
    231         // accumulate data within bin range
    232         for (int j = 0; j < psfTry->sources->n; j++) {
    233             // masked for: bad model fit, outlier in parameters
    234             if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
    235                 continue;
    236 
    237             // skip points with extreme dA values
    238             if (fabs(psfTry->metric->data.F64[j]) > 0.5)
    239                 continue;
    240 
    241             // skip points outside of this bin
    242             if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
    243                 continue;
    244             if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
    245                 continue;
    246 
    247             tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j];
    248             tmp->n ++;
    249         }
    250 
    251         // is this a valid point?
    252         maskB->data.U8[i] = 0;
    253         if (tmp->n < 2) {
    254             maskB->data.U8[i] = 1;
    255             psFree (tmp);
    256             continue;
    257         }
    258 
    259         // dA values are contaminated with low outliers
    260         // measure statistics only on upper 50% of points
    261         // this would be easier if we could sort in reverse:
    262         //
    263         // psVectorSort (tmp, tmp);
    264         // tmp->n = 0.5*tmp->n;
    265         // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    266         // psVectorStats (stats, tmp, NULL, NULL, 0);
    267         // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    268 
    269         psVectorSort (tmp, tmp);
    270         nKeep = 0.5*tmp->n;
    271         nSkip = tmp->n - nKeep;
    272 
    273         psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
    274         for (int j = 0; j < tmp2->n; j++) {
    275             tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
    276         }
    277 
    278         stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    279         psVectorStats (stats, tmp2, NULL, NULL, 0);
    280         psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    281 
    282         daBin->data.F64[i] = stats->sampleMedian;
    283 
    284         psFree (stats);
    285         psFree (tmp);
    286         psFree (tmp2);
    287     }
    288 
    289     // linear fit to rfBin, daBin
    290     psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
    291 
    292     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
    293     poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
    294 
    295     // XXX EAM : replace this when the above version is implemented
    296     // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
    297 
    298     psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    299     psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    300 
    301     stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
    302     stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    303 
    304     psfTry->psf->ApResid = poly->coeff[0];
    305     psfTry->psf->dApResid = stats->clippedStdev;
    306     psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     202    // use 3hi/1lo sigma clipping on the rflux vs metric fit
     203    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     204    stats->min = 1.0;
     205    stats->max = 3.0;
     206    stats->clipIter = 3;
     207
     208    // fit ApResid only to rflux, ignore x,y variations for now
     209    // linear clipped fit of ApResid to rflux
     210    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
     211    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, rflux);
     212    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     213
     214    psfTry->psf->ApTrend->coeff[0][0][0] = poly->coeff[0];
     215    psfTry->psf->ApTrend->coeff[0][0][1] = 0;
     216
     217    psfTry->psf->skyBias  = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     218    psfTry->psf->ApResid  = poly->coeff[0];
     219    psfTry->psf->dApResid = stats->sampleStdev;
    307220
    308221    psFree (rflux);
    309     psFree (rfBin);
    310     psFree (daBin);
    311     psFree (maskB);
    312     psFree (daBinFit);
    313     psFree (daResid);
    314222    psFree (poly);
    315223    psFree (stats);
    … …  
    317225    return true;
    318226}
     227
     228/*
     229  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
     230  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
     231  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
     232*/
     233

  • branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.h

    r5762 r5958  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-12 20:32:44 $
     8 *  @version $Revision: 1.2.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-01-10 04:46:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    109109);
    110110
     111/** pmPSFtryMetric_Alt()
     112 *
     113 * This function is used to measure the PSF model metric for the set of
     114 * results contained in the pmPSFtry structure (alternative implementation).
     115 *
     116 */
     117bool pmPSFtryMetric_Alt(
     118    pmPSFtry *try
     119    ,                      ///< Add comment.
     120    float RADIUS                        ///< Add comment.
     121);
     122
    111123# endif
Note: See TracChangeset for help on using the changeset viewer.