IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28781


Ignore:
Timestamp:
Jul 29, 2010, 2:35:46 PM (16 years ago)
Author:
eugene
Message:

attempting to get a better starting guess; attempting to downweight neighbor detections

Location:
branches/eam_branches/ipp-20100621/psModules/src/objects
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c

    r28656 r28781  
    126126        dPAR[PM_PAR_SKY]  = +1.0;
    127127        dPAR[PM_PAR_I0]   = +f1;
    128         dPAR[PM_PAR_7]    = (z < 0.01) ? -1.5*z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -1.5*z0*f2*log(z);
     128        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
     129        dPAR[PM_PAR_7]   *= 3.0;
    129130
    130131        assert (isfinite(z1));
    … …  
    133134        dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    134135        dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    135         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     136        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
    136137        dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    137138        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    … …  
    224225    psF32     *PAR  = model->params->data.F32;
    225226
     227    // the other parameters depend on the guess for PAR_7
     228    if (!isfinite(PAR[PM_PAR_7])) {
     229        PAR[PM_PAR_7]    = 0.25;
     230    }   
     231    float index = 0.5 / PAR[PM_PAR_7];
     232
     233    // the scale-length is a function of the moments and the index:
     234    // Rmajor_guess = Rmajor_moments * Scale * Zero
     235    float Scale = 0.70 + 0.053 * PAR[PM_PAR_7];
     236    float Zero  = 1.16 - 0.615 * PAR[PM_PAR_7];
     237
    226238    psEllipseMoments emoments;
    227239    emoments.x2 = moments->Mxx;
    … …  
    236248    if (!isfinite(axes.theta)) return false;
    237249
     250    // set a lower limit to avoid absurd solutions..
     251    float Rmajor = PS_MAX(1.0, Scale * axes.major + Zero);
     252    float Rminor = Rmajor * (axes.minor / axes.major);
     253
     254    fprintf (stderr, "guess index: %f : %f, %f -> %f, %f\n", index, axes.major, axes.minor, Rmajor, Rminor);
     255
     256    axes.major = Rmajor;
     257    axes.minor = Rminor;
    238258    psEllipseShape shape = psEllipseAxesToShape (axes);
    239259
    … …  
    242262    if (!isfinite(shape.sxy)) return false;
    243263
    244     // the other parameters depend on the guess for PAR_7
    245     if (!isfinite(PAR[PM_PAR_7])) {
    246         PAR[PM_PAR_7]    = 0.25;
    247         // PAR[PM_PAR_7]    = 0.125;
    248         // PAR[PM_PAR_7]    = 0.5;
    249     }   
    250 
    251     float index = 0.5 / PAR[PM_PAR_7];
    252264    float bn = 1.9992*index - 0.3271;
    253265    // float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
    254266    float Io = exp(0.5*bn);
     267
     268    // XXX do we need this factor of sqrt(2)?
     269    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
     270    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
    255271
    256272    float Sxx = PS_MAX(0.5, shape.sx);
    … …  
    329345    psF64 zn = log(PAR[PM_PAR_I0] / flux);
    330346    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
     347
     348    fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major);
    331349
    332350    psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c

    r28702 r28781  
    256256    return pcm;
    257257}
     258
     259// has the set of fitted terms changed?  has the fitting radius changed?
     260bool pmPCMupdate(pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model) {
     261
     262    bool newWindow = (source->pixels->numRows != pcm->modelFlux->numRows) || (source->pixels->numCols != pcm->modelFlux->numCols);
     263
     264    // re-count the number of unmasked pixels:
     265    if (newWindow) {
     266        for (psS32 i = 0; i < source->pixels->numRows; i++) {
     267            for (psS32 j = 0; j < source->pixels->numCols; j++) {
     268                // XXX are we doing the right thing with the mask?
     269                // skip masked points
     270                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
     271                    continue;
     272                }
     273                // skip zero-variance points
     274                if (source->variance->data.F32[i][j] == 0) {
     275                    continue;
     276                }
     277                // skip nan value points
     278                if (!isfinite(source->pixels->data.F32[i][j])) {
     279                    continue;
     280                }
     281                pcm->nPix ++;
     282            }
     283        }   
     284    }
     285
     286    // if we changed the fit mode, we need to update nDOF
     287    int nParams = 0;
     288    // set parameter mask based on fitting mode
     289    switch (fitOptions->mode) {
     290      case PM_SOURCE_FIT_NORM:
     291        // NORM-only model fits only source normalization (Io)
     292        nParams = 1;
     293        psVectorInit (pcm->constraint->paramMask, 1);
     294        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     295        break;
     296      case PM_SOURCE_FIT_PSF:
     297        // PSF model only fits x,y,Io
     298        nParams = 3;
     299        psVectorInit (pcm->constraint->paramMask, 1);
     300        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     301        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
     302        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
     303        break;
     304      case PM_SOURCE_FIT_EXT:
     305        // EXT model fits all params (except sky)
     306        nParams = model->params->n - 1;
     307        psVectorInit (pcm->constraint->paramMask, 0);
     308        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
     309        break;
     310      case PM_SOURCE_FIT_INDEX:
     311        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
     312        psVectorInit (pcm->constraint->paramMask, 1);
     313        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     314        if (model->params->n == 7) {
     315            nParams = 1;
     316        } else {
     317            nParams = 2;
     318            pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
     319        }
     320        break;
     321      case PM_SOURCE_FIT_NO_INDEX:
     322        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
     323        psVectorInit (pcm->constraint->paramMask, 0);
     324        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
     325        if (model->params->n == 7) {
     326            nParams = model->params->n - 1;
     327        } else {
     328            nParams = model->params->n - 2;
     329            pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
     330        }
     331        break;
     332      default:
     333        psAbort("invalid fitting mode");
     334    }
     335
     336    if (pcm->nPix <  nParams + 1) {
     337        psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
     338        model->flags |= PM_MODEL_STATUS_BADARGS;
     339        return false;
     340    }
     341    pcm->nDOF = pcm->nPix - nParams - 1;
     342
     343    // has the source pixel window changed?
     344    if (newWindow) {
     345
     346        // adjust all supporting images:
     347        pcm->modelFlux = psImageCopy (pcm->modelFlux, source->pixels, PS_TYPE_F32);
     348        for (psS32 n = 0; n < pcm->dmodelsFlux->n; n++) {
     349            if ((pcm->constraint->paramMask != NULL) && (pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     350            pcm->dmodelsFlux->data[n] = psImageCopy (pcm->dmodelsFlux->data[n], source->pixels, PS_TYPE_F32);
     351        }
     352
     353        // adjust images for convolved model and derivative images
     354        pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, source->pixels, PS_TYPE_F32);
     355        for (psS32 n = 0; n < pcm->dmodelsConvFlux->n; n++) {
     356            if ((pcm->constraint->paramMask != NULL) && (pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     357            pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32);
     358        }
     359    }
     360
     361    return true;
     362}

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h

    r28702 r28781  
    6262
    6363pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize);
     64bool pmPCMupdate(pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model);
    6465
    6566psImage *pmPCMdataSaveImage (pmPCMdata *pcm);

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.c

    r28676 r28781  
    966966    bool addNoise = mode & PM_MODEL_OP_NOISE;
    967967
    968     if (source->modelFlux) {
     968    // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
     969    if (!addNoise && source->modelFlux) {
    969970        // add in the pixels from the modelFlux image
    970971        int dX = source->modelFlux->col0 - source->pixels->col0;
    … …  
    987988
    988989        psF32 **target = source->pixels->data.F32;
    989         if (addNoise) {
    990             // when adding noise, we assume the shape and Io have been modified
    991             target = source->variance->data.F32;
    992         }
    993990
    994991        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
    … …  
    10051002            }
    10061003        }
    1007         if (!addNoise) {
    1008             if (add) {
    1009                 source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    1010             } else {
    1011                 source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    1012             }
     1004        if (add) {
     1005            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     1006        } else {
     1007            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    10131008        }
    10141009        return true;
    … …  
    10291024        }
    10301025    }
     1026
     1027    return true;
     1028}
     1029
     1030// should we call pmSourceCacheModel if it does not exist?
     1031bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy)
     1032{
     1033    assert (mode & PM_MODEL_OP_NOISE);
     1034    PS_ASSERT_PTR_NON_NULL(source, false);
     1035    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1036
     1037    if (add) {
     1038        psTrace ("psphot", 3, "adding noise to object at %f,%f\n", source->peak->xf, source->peak->yf);
     1039    } else {
     1040        psTrace ("psphot", 3, "removing noise from object at %f,%f\n", source->peak->xf, source->peak->yf);
     1041    }
     1042
     1043    pmSourceNoiseOpModel (source->modelPSF, source, mode, FACTOR, SIZE, add, maskVal, dx, dy);
     1044
     1045    if (source->modelEXT) {
     1046        pmSourceNoiseOpModel (source->modelEXT, source, mode, FACTOR, SIZE, add, maskVal, dx, dy);
     1047    }
     1048
     1049    return true;
     1050}
     1051
     1052bool pmSourceNoiseOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy)
     1053{
     1054    bool status;
     1055    psEllipseShape oldshape;
     1056    psEllipseShape newshape;
     1057    psEllipseAxes axes;
     1058
     1059    if (add) {
     1060        psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     1061    } else {
     1062        psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     1063    }
     1064
     1065    psF32 *PAR = model->params->data.F32;
     1066
     1067    // save original values
     1068    float oldI0  = PAR[PM_PAR_I0];
     1069    oldshape.sx  = PAR[PM_PAR_SXX];
     1070    oldshape.sy  = PAR[PM_PAR_SYY];
     1071    oldshape.sxy = PAR[PM_PAR_SXY];
     1072
     1073    // XXX can this be done more intelligently?
     1074    if (oldI0 == 0.0) return false;
     1075    if (!isfinite(oldI0)) return false;
     1076
     1077    // increase size and height of source
     1078    axes = psEllipseShapeToAxes (oldshape, 20.0);
     1079    axes.major *= SIZE;
     1080    axes.minor *= SIZE;
     1081    newshape = psEllipseAxesToShape (axes);
     1082    PAR[PM_PAR_I0]  = FACTOR*oldI0;
     1083    PAR[PM_PAR_SXX] = newshape.sx;
     1084    PAR[PM_PAR_SYY] = newshape.sy;
     1085    PAR[PM_PAR_SXY] = newshape.sxy;
     1086
     1087    psImage *target = source->variance;
     1088
     1089    if (add) {
     1090        status = pmModelAddWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
     1091    } else {
     1092        status = pmModelSubWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
     1093    }
     1094
     1095    // restore original values
     1096    PAR[PM_PAR_I0]  = oldI0;
     1097    PAR[PM_PAR_SXX] = oldshape.sx;
     1098    PAR[PM_PAR_SYY] = oldshape.sy;
     1099    PAR[PM_PAR_SXY] = oldshape.sxy;
    10311100
    10321101    return true;

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.h

    r28676 r28781  
    236236bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psImageMaskType maskVal, int dx, int dy);
    237237
     238bool pmSourceNoiseOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
     239bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
     240
    238241bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psImageMaskType maskVal, int dx, int dy);
    239242bool pmSourceCacheModel (pmSource *source, psImageMaskType maskVal);

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c

    r28692 r28781  
    107107bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
    108108
    109     pcm->modelConv->modelGuess(pcm->modelConv, source);
     109    if (!pcm->modelConv->modelGuess(pcm->modelConv, source)) {
     110        return false;
     111    }
     112    return true;
    110113
    111114    // generate copy of the model
    … …  
    114117
    115118    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
     119    // SXX,SYY in model parameters are sqrt(2) * shape.Sxx,Syy
    116120    psEllipseShape psfShape;
    117121    psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
    118122    psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
    119123    psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     124
    120125    psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
     126    if (!isfinite(psfAxes.major)) return false;
     127    if (!isfinite(psfAxes.minor)) return false;
     128    if (!isfinite(psfAxes.theta)) return false;
    121129
    122130    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
    … …  
    125133    extShape.sxy = pcm->modelConv->params->data.F32[PM_PAR_SXY];
    126134    extShape.sy  = pcm->modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     135
    127136    psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0);
     137    if (!isfinite(extAxes.major)) return false;
     138    if (!isfinite(extAxes.minor)) return false;
     139    if (!isfinite(extAxes.theta)) return false;
    128140
    129141    // decrease the initial guess ellipse by psf_minor axis:
    130142    psEllipseAxes extAxesMod;
    131     extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
    132     extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
     143    extAxesMod.major = sqrt (PS_MAX (0.25, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
     144    extAxesMod.minor = sqrt (PS_MAX (0.25, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
    133145    extAxesMod.theta = extAxes.theta;
    134146
    135147    psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod);
     148    if (!isfinite(extShapeMod.sx))  return false;
     149    if (!isfinite(extShapeMod.sy))  return false;
     150    if (!isfinite(extShapeMod.sxy)) return false;
     151
    136152    pcm->modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2;
    137153    pcm->modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy;
Note: See TracChangeset for help on using the changeset viewer.