IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27531


Ignore:
Timestamp:
Mar 30, 2010, 1:29:50 PM (16 years ago)
Author:
eugene
Message:

add diff stats; add flux to CMF_PS1_DV1; output PSF residual mask; more accurate calculation of radius of source model

Location:
trunk/psModules/src
Files:
15 edited
3 copied

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/Makefile.am

    r26162 r27531  
    2222        pmSourceMasks.c \
    2323        pmSourceMoments.c \
     24        pmSourceDiffStats.c \
    2425        pmSourceExtendedPars.c \
    2526        pmSourceUtils.c \
    … …  
    4041        pmSourceIO_CMF_PS1_V1.c \
    4142        pmSourceIO_CMF_PS1_V2.c \
     43        pmSourceIO_CMF_PS1_DV1.c \
    4244        pmSourceIO_MatchedRefs.c \
    4345        pmSourcePlots.c \
    … …  
    7981        pmSource.h \
    8082        pmSourceMasks.h \
     83        pmSourceDiffStats.h \
    8184        pmSourceExtendedPars.h \
    8285        pmSourceUtils.h \

  • trunk/psModules/src/objects/models/pmModel_PGAUSS.c

    r26916 r27531  
    236236psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    237237{
    238     psF64 z, f;
     238    psF64 z;
    239239    int Nstep = 0;
    240240    psEllipseShape shape;
    … …  
    271271    z0 = 0.0;
    272272
    273     // perform a type of bisection to find the value
    274     float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0);
    275     float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0);
    276     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    277         z = 0.5*(z0 + z1);
    278         f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
    279         if (f > limit) {
    280             z0 = z;
    281             f0 = f;
    282         } else {
    283             z1 = z;
    284             f1 = f;
    285         }
    286         Nstep ++;
    287     }
     273    // starting guess:
     274    z = 0.5*(z0 + z1);
     275    float dz = 1.0;
     276
     277    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     278        // use Newton-Raphson to minimize f(z) - limit = 0
     279        float dqdz = (1.0 + z + z*z/2.0);
     280        float q = (dqdz + z*z*z/6.0);
     281
     282        float f = 1.0 / q;
     283        float dfdz = -dqdz * f / q;
     284
     285        dz = (f - limit) / dfdz;
     286
     287        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     288        z -= dz;
     289    }
     290
    288291    psF64 radius = sigma * sqrt (2.0 * z);
    289292

  • trunk/psModules/src/objects/models/pmModel_PS1_V1.c

    r26916 r27531  
    267267psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    268268{
    269     psF64 z, f;
     269    psF64 z;
    270270    int Nstep = 0;
    271271    psEllipseShape shape;
    … …  
    273273    psF32 *PAR = params->data.F32;
    274274
    275     if (flux <= 0) {
    276         return 1.0;
    277     }
    278     if (PAR[PM_PAR_I0] <= 0) {
    279         return 1.0;
    280     }
    281     if (flux >= PAR[PM_PAR_I0]) {
    282         return 1.0;
    283     }
    284     if (PAR[PM_PAR_7] == 0.0) {
    285         return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    286     }
     275    if (flux <= 0) return 1.0;
     276    if (PAR[PM_PAR_I0] <= 0) return 1.0;
     277    if (flux >= PAR[PM_PAR_I0]) return 1.0;
     278    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    287279
    288280    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    … …  
    305297    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
    306298
    307     // perform a type of bisection to find the value
    308     float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, ALPHA));
    309     float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, ALPHA));
    310     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    311         z = 0.5*(z0 + z1);
    312         f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    313         if (f > limit) {
    314             z0 = z;
    315             f0 = f;
    316         } else {
    317             z1 = z;
    318             f1 = f;
    319         }
    320         Nstep ++;
     299    // starting guess:
     300    z = 0.5*(z0 + z1);
     301    float dz = 1.0;
     302
     303    // use Newton-Raphson to minimize f(z) - limit = 0
     304    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     305        float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     306        float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
     307
     308        float f = 1.0 / q;
     309        float dfdz = -dqdz * f / q;
     310
     311        dz = (f - limit) / dfdz;
     312
     313        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     314        z -= dz;
    321315    }
    322316    psF64 radius = sigma * sqrt (2.0 * z);
    … …  
    350344    // convert to shape terms (SXX,SYY,SXY)
    351345    if (!pmPSF_FitToModel (out, 0.1)) {
    352         // psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    353346        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    354347        return false;
    … …  
    475468}
    476469
    477 
    478470# undef PM_MODEL_FUNC
    479471# undef PM_MODEL_FLUX

  • trunk/psModules/src/objects/models/pmModel_QGAUSS.c

    r26916 r27531  
    3939# define PM_MODEL_SET_LIMITS      pmModelSetLimits_QGAUSS
    4040
     41# define ALPHA   2.250
     42# define ALPHA_M 1.250
     43
    4144// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    4245// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
    … …  
    4447
    4548// Lax parameter limits
    46 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.1 };
     49static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
    4750static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
    4851
    4952// Moderate parameter limits
    50 static float *paramsMinModerate = paramsMinLax;
    51 static float *paramsMaxModerate = paramsMaxLax;
     53// Tolerate a small divot (k < 0)
     54static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
     55static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
    5256
    5357// Strict parameter limits
    54 static float *paramsMinStrict = paramsMinLax;
    55 static float *paramsMaxStrict = paramsMaxLax;
     58// k = PAR_7 < 0 is very undesirable (big divot in the middle)
     59static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
     60static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
    5661
    5762// Parameter limits to use
    5863static float *paramsMinUse = paramsMinLax;
    5964static float *paramsMaxUse = paramsMaxLax;
    60 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
     65static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
    6166
    6267static bool limitsApply = true;         // Apply limits?
    … …  
    8186    assert (z >= 0);
    8287
    83     psF32 zp = pow(z,1.25);
     88    psF32 zp = pow(z,ALPHA_M);
    8489    psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp);
    8590
    … …  
    9297        // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r
    9398        psF32 t = r1*r;
    94         psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp);
     99        psF32 q = t*(PAR[PM_PAR_7] + ALPHA*zp);
    95100
    96101        dPAR[PM_PAR_SKY]  = +1.0;
    … …  
    246251    float f1, f2;
    247252    for (z = DZ; z < 50; z += DZ) {
    248         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     253        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    249254        z += DZ;
    250         f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     255        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    251256        norm += f0 + 4*f1 + f2;
    252257        f0 = f2;
    … …  
    263268psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    264269{
    265     psF64 z, f;
     270    psF64 z;
    266271    int Nstep = 0;
    267272    psEllipseShape shape;
    … …  
    269274    psF32 *PAR = params->data.F32;
    270275
    271     if (flux <= 0)
    272         return (1.0);
    273     if (PAR[PM_PAR_I0] <= 0)
    274         return (1.0);
    275     if (flux >= PAR[PM_PAR_I0])
    276         return (1.0);
     276    if (flux <= 0) return 1.0;
     277    if (PAR[PM_PAR_I0] <= 0) return 1.0;
     278    if (flux >= PAR[PM_PAR_I0]) return 1.0;
     279    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    277280
    278281    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    … …  
    290293
    291294    // choose a z value guaranteed to be beyond our limit
    292     float z0 = pow((1.0 / limit), (1.0 / 2.25));
    293     psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]);
    294     float z1 = (1.0 / limit) / PAR[PM_PAR_7];
     295    float z0 = 0.0;
     296    float z1 = pow((1.0 / limit), (1.0 / ALPHA));
    295297    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
    296     z1 = PS_MAX (z0, z1);
    297     z0 = 0.0;
    298 
    299     // perform a type of bisection to find the value
    300     float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25));
    301     float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25));
    302     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    303         z = 0.5*(z0 + z1);
    304         f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    305         if (f > limit) {
    306             z0 = z;
    307             f0 = f;
    308         } else {
    309             z1 = z;
    310             f1 = f;
    311         }
    312         Nstep ++;
     298    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
     299
     300    // starting guess:
     301    z = 0.5*(z0 + z1);
     302    float dz = 1.0;
     303
     304    // use Newton-Raphson to minimize f(z) - limit = 0
     305    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     306        float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     307        float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
     308
     309        float f = 1.0 / q;
     310        float dfdz = -dqdz * f / q;
     311
     312        dz = (f - limit) / dfdz;
     313
     314        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     315        z -= dz;
    313316    }
    314317    psF64 radius = sigma * sqrt (2.0 * z);
    … …  
    475478# undef PM_MODEL_FIT_STATUS
    476479# undef PM_MODEL_SET_LIMITS
     480# undef ALPHA
     481# undef ALPHA_M

  • trunk/psModules/src/objects/models/pmModel_RGAUSS.c

    r26916 r27531  
    257257psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    258258{
    259     psF64 z, f;
     259    psF64 z;
    260260    int Nstep = 0;
    261261    psEllipseShape shape;
    … …  
    291291    z0 = 0.0;
    292292
    293     // perform a type of bisection to find the value
    294     float f0 = 1.0 / (1 + z0 + pow(z0, PAR[PM_PAR_7]));
    295     float f1 = 1.0 / (1 + z1 + pow(z1, PAR[PM_PAR_7]));
    296     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    297         z = 0.5*(z0 + z1);
    298         f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
    299         if (f > limit) {
    300             z0 = z;
    301             f0 = f;
    302         } else {
    303             z1 = z;
    304             f1 = f;
    305         }
    306         Nstep ++;
    307     }
     293    // starting guess:
     294    z = 0.5*(z0 + z1);
     295    float dz = 1.0;
     296
     297    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     298        // use Newton-Raphson to minimize f(z) - limit = 0
     299        float q = (1 + z + pow(z,PAR[PM_PAR_7]));
     300        float dqdz = (1.0 + PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0));
     301
     302        float f = 1.0 / q;
     303        float dfdz = -dqdz * f / q;
     304
     305        dz = (f - limit) / dfdz;
     306
     307        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     308        z -= dz;
     309    }
     310
    308311    psF64 radius = sigma * sqrt (2.0 * z);
    309312

  • trunk/psModules/src/objects/pmGrowthCurveGenerate.c

    r25754 r27531  
    157157
    158158    // measure the fitMag for this model
    159     pmSourcePhotometryModel (&fitMag, model);
     159    pmSourcePhotometryModel (&fitMag, NULL, model);
    160160    growth->fitMag = fitMag;
    161161

  • trunk/psModules/src/objects/pmPSF_IO.c

    r27178 r27531  
    491491        psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
    492492
    493         // write the residuals as three planes of the image
    494         // this call creates an extension with NAXIS3 = 3
     493        // write the residuals as planes of the image
     494        psArray *images = psArrayAllocEmpty (1);
     495        psArrayAdd (images, 1, psf->residuals->Ro);  // z = 0 is Ro
     496
    495497        if (psf->residuals->Rx) {
    496             // this call creates an extension with NAXIS3 = 3
    497             psArray *images = psArrayAllocEmpty (3);
    498             psArrayAdd (images, 1, psf->residuals->Ro);
    499498            psArrayAdd (images, 1, psf->residuals->Rx);
    500499            psArrayAdd (images, 1, psf->residuals->Ry);
    501 
    502             if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
    503                 psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
    504                 psFree(images);
    505                 psFree(residName);
    506                 psFree(header);
    507                 return false;
    508             }
    509             psFree (images);
    510         } else {
    511             // this call creates an extension with NAXIS3 = 1
    512             if (!psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName)) {
    513                 psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
    514                 psFree(residName);
    515                 psFree(header);
    516                 return false;
    517             }
    518         }
     500        }
     501
     502        // note that all N plane are implicitly of the same type, so we convert the mask
     503        if (psf->residuals->mask) {
     504            psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
     505            psArrayAdd (images, 1, mask);
     506            psFree (mask);
     507        }
     508
     509        // psFitsWriteImageCube (file->fits, header, images, residName);
     510        // psFree (images);
     511
     512        if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
     513            psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
     514            psFree(images);
     515            psFree(residName);
     516            psFree(header);
     517            return false;
     518        }
     519        psFree (images);
    519520        psFree (residName);
    520521        psFree (header);
    … …  
    10171018            return false;
    10181019        }
    1019         if (Nz > 1) {
    1020             assert (Nz == 3);
     1020
     1021        // note that all N plane are implicitly of the same type, so we convert the mask
     1022        psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
     1023        psImageInit (psf->residuals->mask, 0);
     1024        psImageInit (psf->residuals->Rx, 0.0);
     1025        psImageInit (psf->residuals->Ry, 0.0);
     1026        switch (Nz) {
     1027          case 1: // Ro only
     1028            break;
     1029          case 2: // Ro and mask
     1030            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 1)) {
     1031                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
     1032                return false;
     1033            }
     1034            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
     1035            break;
     1036          case 3: // Ro, Rx and Ry, no mask
    10211037            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
    10221038                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
    … …  
    10271043                return false;
    10281044            }
    1029         }
    1030         // XXX notice that we are not saving the resid->mask
     1045            break;
     1046          case 4: // Ro, Rx, Ry, and mask:
     1047            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
     1048                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
     1049                return false;
     1050            }
     1051            if (!psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2)) {
     1052                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
     1053                return false;
     1054            }
     1055            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 3)) {
     1056                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
     1057                return false;
     1058            }
     1059            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
     1060            break;
     1061        }
     1062        psFree (mask);
    10311063    }
    10321064

  • trunk/psModules/src/objects/pmSource.c

    r26893 r27531  
    4646    psFree(tmp->maskView);
    4747    psFree(tmp->modelFlux);
    48     psFree(tmp->psfFlux);
     48    psFree(tmp->psfImage);
    4949    psFree(tmp->moments);
    5050    psFree(tmp->modelPSF);
    … …  
    5252    psFree(tmp->modelFits);
    5353    psFree(tmp->extpars);
     54    psFree(tmp->moments);
     55    psFree(tmp->diffStats);
    5456    psFree(tmp->blends);
    5557    psTrace("psModules.objects", 10, "---- end ----\n");
    … …  
    6870    psFree (source->maskView);
    6971    psFree (source->modelFlux);
    70     psFree (source->psfFlux);
     72    psFree (source->psfImage);
    7173
    7274    source->pixels = NULL;
    … …  
    7577    source->maskView = NULL;
    7678    source->modelFlux = NULL;
    77     source->psfFlux = NULL;
     79    source->psfImage = NULL;
    7880    return;
    7981}
    … …  
    103105    source->maskView = NULL;
    104106    source->modelFlux = NULL;
    105     source->psfFlux = NULL;
     107    source->psfImage = NULL;
    106108    source->moments = NULL;
    107109    source->blends = NULL;
    … …  
    113115    source->tmpFlags = 0;
    114116    source->extpars = NULL;
     117    source->diffStats = NULL;
     118
    115119    source->region = psRegionSet(NAN, NAN, NAN, NAN);
    116120    psMemSetDeallocator(source, (psFreeFunc) sourceFree);
    117121
    118122    // default values are NAN
    119     source->psfMag = NAN;
     123    source->psfMag     = NAN;
     124    source->psfFlux    = NAN;
     125    source->psfFluxErr = NAN;
    120126    source->extMag = NAN;
    121127    source->errMag = NAN;
    … …  
    259265        mySource->modelFlux = NULL;
    260266
    261         // drop the old psfFlux pixels and force the user to re-create
    262         psFree (mySource->psfFlux);
    263         mySource->psfFlux = NULL;
     267        // drop the old psfImage pixels and force the user to re-create
     268        psFree (mySource->psfImage);
     269        mySource->psfImage = NULL;
    264270    }
    265271    return extend;
    … …  
    873879
    874880    // if we already have a cached image, re-use that memory
    875     source->psfFlux = psImageCopy (source->psfFlux, source->pixels, PS_TYPE_F32);
    876     psImageInit (source->psfFlux, 0.0);
     881    source->psfImage = psImageCopy (source->psfImage, source->pixels, PS_TYPE_F32);
     882    psImageInit (source->psfImage, 0.0);
    877883
    878884    // in some places (psphotEnsemble), we need a normalized version
    879885    // in others, we just want the model.  which is more commonly used?
    880     // psfFlux always has unity normalization (I0 = 1.0)
    881     pmModelAdd (source->psfFlux, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     886    // psfImage always has unity normalization (I0 = 1.0)
     887    pmModelAdd (source->psfImage, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
    882888    return true;
    883889}

  • trunk/psModules/src/objects/pmSource.h

    r26893 r27531  
    1616#include "pmMoments.h"
    1717#include "pmSourceExtendedPars.h"
     18#include "pmSourceDiffStats.h"
    1819
    1920/// @addtogroup Objects Object Detection / Analysis Functions
    … …  
    3839
    3940typedef enum {
    40     PM_SOURCE_TMPF_MODEL_GUESS   = 0x0001,
    41     PM_SOURCE_TMPF_SUBTRACTED    = 0x0002,
    42     PM_SOURCE_TMPF_SIZE_MEASURED = 0x0004,
     41    PM_SOURCE_TMPF_MODEL_GUESS       = 0x0001,
     42    PM_SOURCE_TMPF_SUBTRACTED        = 0x0002,
     43    PM_SOURCE_TMPF_SIZE_MEASURED     = 0x0004,
     44    PM_SOURCE_TMPF_SIZE_CR_CANDIDATE = 0x0008,
    4345} pmSourceTmpF;
    4446
    … …  
    6466    psImage *maskView;                  ///< view into global image mask for this object region
    6567    psImage *modelFlux;                 ///< cached copy of the best model for this source
    66     psImage *psfFlux;                   ///< cached copy of the psf model for this source
     68    psImage *psfImage;                   ///< cached copy of the psf model for this source
    6769    pmMoments *moments;                 ///< Basic moments measured for the object.
    6870    pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
    … …  
    7476    psArray *blends;                    ///< collection of sources thought to be confused with object
    7577    float psfMag;                       ///< calculated from flux in modelPSF
     78    float psfFlux;                      ///< calculated from flux in modelPSF
     79    float psfFluxErr;                   ///< calculated from flux in modelPSF
    7680    float extMag;                       ///< calculated from flux in modelEXT
    7781    float errMag;                       ///< error in psfMag OR extMag (depending on type)
    … …  
    8589    psRegion region;                    ///< area on image covered by selected pixels
    8690    pmSourceExtendedPars *extpars;      ///< extended source parameters
     91    pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
    8792};
    8893

  • trunk/psModules/src/objects/pmSourceIO.c

    r26893 r27531  
    539539                status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
    540540            }
     541            if (!strcmp (exttype, "PS1_DV1")) {
     542                status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname);
     543            }
    541544
    542545            if (xsrcname) {
    … …  
    553556                    status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
    554557                }
     558                if (!strcmp (exttype, "PS1_DV1")) {
     559                    status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, sources, xsrcname, recipe);
     560                }
    555561            }
    556562            if (xfitname) {
    … …  
    567573                    status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
    568574                }
     575                if (!strcmp (exttype, "PS1_DV1")) {
     576                    status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, sources, xfitname);
     577                }
    569578            }
    570579            psFree (outhead);
    … …  
    10151024                sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
    10161025            }
     1026            if (!strcmp (exttype, "PS1_DV1")) {
     1027                sources = pmSourcesRead_CMF_PS1_DV1 (file->fits, hdu->header);
     1028            }
    10171029
    10181030            if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {

  • trunk/psModules/src/objects/pmSourceIO.h

    r24694 r27531  
    4343bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname);
    4444
     45bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
     46bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
     47bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname);
     48
    4549bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    4650
    … …  
    5357psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
    5458psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
     59psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
    5560
    5661bool pmSourcesWritePSFs (psArray *sources, char *filename);

  • trunk/psModules/src/objects/pmSourceMoments.c

    r26893 r27531  
    163163    }
    164164
    165     // if we have less than (1/2) of the possible pixels, force a retry
     165    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
     166    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
     167
    166168    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    167     if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    168         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
     169    if ((numPixels < minPixels) || (Sum <= 0)) {
     170        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
    169171        return (false);
    170172    }

  • trunk/psModules/src/objects/pmSourcePhotometry.c

    r25980 r27531  
    109109        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
    110110        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
    111         source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
     111        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
     112        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
     113        source->psfMag = -2.5*log10(source->psfFlux);
    112114    } else {
    113         status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
     115        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
     116        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
    114117    }
    115118
    … …  
    119122        for (int i = 0; i < source->modelFits->n; i++) {
    120123            pmModel *model = source->modelFits->data[i];
    121             status = pmSourcePhotometryModel (&model->mag, model);
     124            status = pmSourcePhotometryModel (&model->mag, NULL, model);
    122125            if (model == source->modelEXT) foundEXT = true;
    123126        }
    … …  
    125128            source->extMag = source->modelEXT->mag;
    126129        } else {
    127             status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
     130            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
    128131        }
    129132    } else {
    130133        if (source->modelEXT) {
    131             status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
     134            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
    132135        }
    133136    }
    … …  
    143146    if (mode & PM_SOURCE_PHOT_WEIGHT) {
    144147        pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
     148    }
     149
     150    // measure the contribution of included pixels
     151    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
     152        pmSourceMeasureDiffStats (source, maskVal);
    145153    }
    146154
    … …  
    217225
    218226// return source model magnitude
    219 bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
    220 {
    221     PS_ASSERT_PTR_NON_NULL(fitMag, false);
    222     if (model == NULL) {
    223         return false;
    224     }
    225 
    226     float fitSum = 0;
    227     *fitMag = NAN;
     227bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model)
     228{
     229    psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)");
     230    if (model == NULL) return false;
     231
     232    float mag  = NAN;
     233    float flux = NAN;
    228234
    229235    // measure fitMag
    230     fitSum = model->modelFlux (model->params);
    231     if (fitSum <= 0)
    232         return false;
    233     if (!isfinite(fitSum))
    234         return false;
    235     *fitMag = -2.5*log10(fitSum);
     236    flux = model->modelFlux (model->params);
     237    if (flux > 0) {
     238        mag = -2.5*log10(flux);
     239    }
     240    if (fitMag) {
     241        *fitMag = mag;
     242    }
     243    if (fitFlux) {
     244        *fitFlux = flux;
     245    }
     246
     247    if (flux <= 0) return false;
     248    if (!isfinite(flux)) return false;
    236249
    237250    return (true);
    … …  
    356369
    357370    *pixWeight = validSum / modelSum;
     371    return (true);
     372}
     373
     374# define FLUX_LIMIT 3.0
     375
     376// return source aperture magnitude
     377bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
     378{
     379    PS_ASSERT_PTR_NON_NULL(source, false);
     380
     381    if (source->diffStats == NULL) {
     382        source->diffStats = pmSourceDiffStatsAlloc();
     383    }
     384
     385    float fGood = 0.0;
     386    float fBad  = 0.0;
     387    int   nGood = 0;
     388    int   nMask = 0;
     389    int   nBad  = 0;
     390   
     391    psImage *flux     = source->pixels;
     392    psImage *variance = source->variance;
     393    psImage *mask     = source->maskObj;
     394
     395    for (int iy = 0; iy < flux->numRows; iy++) {
     396        for (int ix = 0; ix < flux->numCols; ix++) {
     397            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
     398                nMask ++;
     399                continue;
     400            }
     401
     402            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     403
     404            if (SN > +FLUX_LIMIT) {
     405                nGood ++;
     406                fGood += fabs(flux->data.F32[iy][ix]);
     407            }
     408
     409            if (SN < -FLUX_LIMIT) {
     410                nBad ++;
     411                fBad += fabs(flux->data.F32[iy][ix]);
     412            }
     413        }
     414    }
     415
     416    source->diffStats->nGood      = nGood;
     417    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);         
     418    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);         
     419    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);         
     420    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
     421
    358422    return (true);
    359423}

  • trunk/psModules/src/objects/pmSourcePhotometry.h

    r25980 r27531  
    2929
    3030typedef enum {
    31     PM_SOURCE_PHOT_NONE   = 0x0000,
    32     PM_SOURCE_PHOT_GROWTH = 0x0001,
    33     PM_SOURCE_PHOT_APCORR = 0x0002,
    34     PM_SOURCE_PHOT_WEIGHT = 0x0004,
    35     PM_SOURCE_PHOT_INTERP = 0x0008,
     31    PM_SOURCE_PHOT_NONE      = 0x0000,
     32    PM_SOURCE_PHOT_GROWTH    = 0x0001,
     33    PM_SOURCE_PHOT_APCORR    = 0x0002,
     34    PM_SOURCE_PHOT_WEIGHT    = 0x0004,
     35    PM_SOURCE_PHOT_INTERP    = 0x0008,
     36    PM_SOURCE_PHOT_DIFFSTATS = 0x0010,
    3637} pmSourcePhotometryMode;
    3738
    3839bool pmSourcePhotometryModel(
    3940    float *fitMag,                      ///< integrated fit magnitude
     41    float *fitFlux,                     ///< integrated fit magnitude
    4042    pmModel *model                      ///< model used for photometry
    4143);
    … …  
    5456bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor);
    5557
     58bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal);
    5659
    5760double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);

  • trunk/psModules/src/psmodules.h

    r26893 r27531  
    115115#include <pmDetections.h>
    116116#include <pmMoments.h>
     117#include <pmSourceExtendedPars.h>
     118#include <pmSourceDiffStats.h>
    117119#include <pmResiduals.h>
    118120#include <pmGrowthCurve.h>
Note: See TracChangeset for help on using the changeset viewer.