IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 10, 2013, 2:55:11 PM (13 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20130904

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psphot

  • trunk/psphot/src

  • trunk/psphot/src/psphotSourceFits.c

    r36107 r36375  
    1919static int NfitIterPCM = 0;
    2020static int NfitPixPCM = 0;
     21
     22bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize);
     23bool psphotPCMfitRetry (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, float psfSize);
    2124
    2225bool psphotFitInit (int nThreads) {
    … …  
    579582    maskVal |= markVal;
    580583
    581     // allocate the model
     584    // allocate the model (this can only fail on a config error)
    582585    pmModel *model = pmModelAlloc(modelType);
    583     if (!model) {
    584         return NULL;
    585     }
    586 
    587     float t1, t2, t4, t5;
    588     t1 = t2 = t4 = t5 = 0.0;
     586    psAssert (model, "invalid extended model name");
     587
     588    float t1, t2, t3, t4, t5;
     589    t1 = t2 = t3 = t4 = t5 = 0.0;
    589590    if (TIMING) { psTimerStart ("psphotFitPCM"); }
    590591
    … …  
    619620        if (!psphotSersicModelGuessPCM (pcm, source, maskVal, psfSize)) {
    620621            psFree (pcm);
    621             model->flags |= PM_MODEL_STATUS_BADARGS;
     622            model->flags |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
    622623            return model;
    623624        }
    … …  
    628629        if (!psphotFitSersicShapeAndIndexGridAuto (pcm, readout, source, &options, maskVal, markVal, psfSize)) {
    629630            psFree (pcm);
    630             model->flags |= PM_MODEL_STATUS_BADARGS;
     631            model->flags |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
    631632            psError(PS_ERR_UNKNOWN, true, "Failed to find a index & shape");
    632633            psErrorClear (); // clear the polynomial error
    … …  
    637638        if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
    638639            psFree (pcm);
    639             model->flags |= PM_MODEL_STATUS_BADARGS;
     640            model->flags |= PM_MODEL_PCM_FAIL_GUESS;
    640641            return model;
    641642        }
    642643    }
    643644
    644     if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
     645    if (TIMING) { t3 = psTimerMark ("psphotFitPCM"); }
    645646
    646647    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    … …  
    648649    NfitIterPCM += pcm->modelConv->nIter;
    649650    NfitPixPCM += pcm->modelConv->nDOF;
     651    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
     652
     653    // XXX we might make this more efficient by setting NITER to be fairly small.  if we hit the iteration
     654    // limit, then we could do a small grid search on the size and try again from the best fit
     655
     656    if (options.isInteractive) psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
     657    // if (pcm->modelConv->nIter == fitOptions->nIter) {
     658    //  psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
     659    // }
    650660    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
    651661
    652662    if (TIMING) {
    653663        int nPixBig = source->pixels->numCols * source->pixels->numRows;
    654         fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
     664        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t3: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t3, t4, t5);
    655665    }
    656666    if (EXTRA_VERBOSE && !TIMING) {
    … …  
    778788}
    779789
    780 // float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
    781 float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
    782 # define N_INDEX_GUESS_INV 7
    783 
    784 // float reffGuess[] = {3.0, 10.0, 20.0, 30.0, 40.0};
    785 float reffGuess[] = {0.5, 0.75, 1.0, 1.4, 2.0};
     790// 0.5 / n for (1.0, 1.25, 1.66, 2.0, 3.33, 4.0)
     791// float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
     792
     793// 0.5 / n for (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0)
     794float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 0.166, 0.125, 0.10, 0.0833};
     795float indexGuessR1q[] = {1.06, 1.19, 1.335, 1.48, 1.840, 2.290, 2.84, 3.5300};
     796# define N_INDEX_GUESS_INV 8
     797
     798// we are going to guess in fractions about the R1-based guess
     799float reffGuess[] = {0.8, 0.9, 1.0, 1.12, 1.25};
    786800# define N_REFF_GUESS 5
    787801
    … …  
    791805
    792806    // we get a reasonable guess from:
    793     // * Reff = Kron R1
     807    // * Reff = Kron R1 / Q(index) -- Q comes from Graham & Driver
    794808    // * Rmajor / Rminor & Theta from moments
    795809    // * Io from total Kron flux
    … …  
    813827    psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0);
    814828
    815     if (0) {
    816         psEllipseAxes guessAxes;
    817         guessAxes.major = source->moments->Mrf;
    818         guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
    819         guessAxes.theta = momentAxes.theta;
    820 
    821         if (!isfinite(guessAxes.major)) return false;
    822         if (!isfinite(guessAxes.minor)) return false;
    823         if (!isfinite(guessAxes.theta)) return false;
    824 
    825         // convert the major,minor,theta to shape parameters for an Reff-like model
    826         pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
    827     }
    828 
    829829    // set the model position
    830830    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
    … …  
    846846    // loop over index guesses and find the best fit
    847847    for (int j = 0; j < N_REFF_GUESS; j++) {
    848 
    849         psEllipseAxes guessAxes;
    850         guessAxes.major = reffGuess[j] * source->moments->Mrf;
    851         guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
    852         guessAxes.theta = momentAxes.theta;
    853 
    854         if (!isfinite(guessAxes.major)) return false;
    855         if (!isfinite(guessAxes.minor)) return false;
    856         if (!isfinite(guessAxes.theta)) return false;
    857 
    858         // convert the major,minor,theta to shape parameters for an Reff-like model
    859         pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
    860 
    861848        for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
    862849            PAR[PM_PAR_7] = indexGuessInv[i];
    863850
     851            psEllipseAxes guessAxes;
     852            guessAxes.major = reffGuess[j] * source->moments->Mrf / indexGuessR1q[i];
     853            guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
     854            guessAxes.theta = momentAxes.theta;
     855
     856            if (!isfinite(guessAxes.major)) return false;
     857            if (!isfinite(guessAxes.minor)) return false;
     858            if (!isfinite(guessAxes.theta)) return false;
     859
     860            // convert the major,minor,theta to shape parameters for an Reff-like model
     861            pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
     862
    864863            // generated the modelFlux
    865             pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
     864            // XXX note that this does not add sky to model
     865            pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
    866866       
    867867            float YY = 0.0;
    … …  
    901901                iMin = Io;
    902902                sMin = indexGuessInv[i];
    903                 rMin = reffGuess[j];
     903                rMin = reffGuess[j] / indexGuessR1q[i];
    904904            }
    905             // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
    906             // fprintf (stderr, "\n");
     905            if (EXTRA_VERBOSE) {
     906                fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
     907                fprintf (stderr, "\n");
     908            }
    907909        }
    908910    }
    … …  
    10071009    psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
    10081010
    1009     PAR[PM_PAR_7] = indexGuessInv[nStart];
     1011    float Sm = NAN, Sp = NAN, So = NAN;
     1012    if (nStart == 0) {
     1013        Sm = indexGuessInv[nStart];
     1014        So = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
     1015        Sp = indexGuessInv[nStart + 1];
     1016    } else if (nStart == N_INDEX_GUESS_INV - 1) {
     1017        Sp = indexGuessInv[nStart];
     1018        So = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
     1019        Sm = indexGuessInv[nStart - 1];
     1020    } else {
     1021        Sm = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
     1022        So = indexGuessInv[nStart];
     1023        Sp = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
     1024    }
     1025   
     1026    PAR[PM_PAR_7] = Sm;
    10101027    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    10111028    if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
    … …  
    10131030    psVectorAppend (chi2, model->chisqNorm);
    10141031
    1015     PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.1;
     1032    PAR[PM_PAR_7] = So;
    10161033    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    10171034    if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
    … …  
    10191036    psVectorAppend (chi2, model->chisqNorm);
    10201037
    1021     PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55;
     1038    PAR[PM_PAR_7] = Sp;
    10221039    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    10231040    if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
    10241041    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
    10251042    psVectorAppend (chi2, model->chisqNorm);
    1026 
    1027     if (chi2->data.F32[1] < chi2->data.F32[2]) {
    1028       if (nStart == N_INDEX_GUESS_INV - 1) {
    1029         PAR[PM_PAR_7] = 0.11;
    1030       } else {
    1031         PAR[PM_PAR_7] = indexGuessInv[nStart + 1];
    1032       }
    1033     } else {
    1034       if (nStart == 0) {
    1035         PAR[PM_PAR_7] = 0.52;
    1036       } else {
    1037         PAR[PM_PAR_7] = indexGuessInv[nStart - 1];
    1038       }
    1039     }
    10401043
    10411044    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
    … …  
    10511054    float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0;
    10521055
    1053     // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev)
    1054     Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1);
     1056    // constrain Smin to be in a valid range: allow the fitted range to go a bit beyond the 3 trial points, but no further
     1057    float Smx = Sm - 0.25*(So - Sm);
     1058    float Spx = Sp + 0.25*(Sp - So);
     1059    Smin = PS_MAX(PS_MIN(Smin, Smx), Spx);
    10551060    PAR[PM_PAR_7] = Smin;
     1061
     1062    // XXX I could set the error on PAR_7 here if I knew how to roughly convert these chisq values to true chisq values
    10561063
    10571064    // return to the original fitting mode (fitOptions)
    … …  
    12461253}
    12471254
    1248  
     1255// # define N_REFF_CHECK 11
     1256// float drefCheck[] = {-0.02, -0.04, -0.06, 0.0, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25};
     1257
     1258// we have an initial fit, check to see if the current size is besst
     1259bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) {
     1260
     1261    // PAR is already at my current best guess
     1262    psF32 *PAR = pcm->modelConv->params->data.F32;
     1263
     1264    // store best guess as a shape
     1265    psEllipseAxes centerAxes;
     1266    pmModelParamsToAxes (&centerAxes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     1267
     1268    float xMin = NAN;
     1269    float iMin = NAN;
     1270    float rMin = NAN;
     1271
     1272    // loop over Reff, keeping the ARatio and Theta constant
     1273    for (int j = -4; j <= 4; j++) {
     1274
     1275        float dref = j * 0.01;
     1276
     1277        psEllipseAxes guessAxes;
     1278        guessAxes.major = centerAxes.major + dref;
     1279        guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
     1280        guessAxes.theta = centerAxes.theta;
     1281
     1282        if (!isfinite(guessAxes.major)) return false;
     1283        if (!isfinite(guessAxes.minor)) return false;
     1284        if (!isfinite(guessAxes.theta)) return false;
     1285
     1286        // convert the major,minor,theta to shape parameters for an Reff-like model
     1287        pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
     1288
     1289        // generated the modelFlux
     1290            // XXX note that this does not add sky to model
     1291        pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
     1292       
     1293        float YY = 0.0;
     1294        float YM = 0.0;
     1295        float MM = 0.0;
     1296        bool usePoisson = false;
     1297
     1298        for (int iy = 0; iy < source->pixels->numRows; iy++) {
     1299            for (int ix = 0; ix < source->pixels->numCols; ix++) {
     1300                // skip masked points
     1301                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     1302                    continue;
     1303                }
     1304                // skip zero-variance points
     1305                if (source->variance->data.F32[iy][ix] == 0) {
     1306                    continue;
     1307                }
     1308                // skip nan value points
     1309                if (!isfinite(source->pixels->data.F32[iy][ix])) {
     1310                    continue;
     1311                }
     1312
     1313                float fy = source->pixels->data.F32[iy][ix];
     1314                float fm = source->modelFlux->data.F32[iy][ix];
     1315                float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
     1316
     1317                YY += PS_SQR(fy) * wt;
     1318                YM += fm * fy * wt;
     1319                MM += PS_SQR(fm) * wt;
     1320            }
     1321        }
     1322
     1323        float Io = YM / MM;
     1324        float Chisq = YY - 2 * Io * YM + Io * Io * MM;
     1325        if (isnan(xMin) || (Chisq < xMin)) {
     1326            xMin = Chisq;
     1327            iMin = Io;
     1328            rMin = dref;
     1329        }
     1330        // fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
     1331    }
     1332
     1333    psEllipseAxes guessAxes;
     1334    guessAxes.major = centerAxes.major + rMin;
     1335    guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
     1336    guessAxes.theta = centerAxes.theta;
     1337
     1338    if (!isfinite(guessAxes.major)) return false;
     1339    if (!isfinite(guessAxes.minor)) return false;
     1340    if (!isfinite(guessAxes.theta)) return false;
     1341
     1342    // convert the major,minor,theta to shape parameters for an Reff-like model
     1343    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
     1344    PAR[PM_PAR_I0] = iMin;
     1345
     1346    return true;
     1347}
     1348
     1349// we have an initial fit, check to see if the current size is besst
     1350bool psphotPCMfitRetry (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, float psfSize) {
     1351
     1352    // PAR is already at my current best guess
     1353    psF32 *PAR = pcm->modelConv->params->data.F32;
     1354
     1355    // store best guess as a shape
     1356    psEllipseAxes centerAxes;
     1357    pmModelParamsToAxes (&centerAxes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     1358
     1359    // retry with axes smaller by 1 pixel
     1360    psEllipseAxes guessAxes;
     1361    guessAxes.major = centerAxes.major - 0.08;
     1362    guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
     1363    guessAxes.theta = centerAxes.theta;
     1364
     1365    if (!isfinite(guessAxes.major)) return false;
     1366    if (!isfinite(guessAxes.minor)) return false;
     1367    if (!isfinite(guessAxes.theta)) return false;
     1368
     1369    // convert the major,minor,theta to shape parameters for an Reff-like model
     1370    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
     1371
     1372    // generated the modelFlux
     1373            // XXX note that this does not add sky to model
     1374    pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
     1375       
     1376    float YY = 0.0;
     1377    float YM = 0.0;
     1378    float MM = 0.0;
     1379    bool usePoisson = false;
     1380
     1381    for (int iy = 0; iy < source->pixels->numRows; iy++) {
     1382        for (int ix = 0; ix < source->pixels->numCols; ix++) {
     1383            // skip masked points
     1384            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     1385                continue;
     1386            }
     1387            // skip zero-variance points
     1388            if (source->variance->data.F32[iy][ix] == 0) {
     1389                continue;
     1390            }
     1391            // skip nan value points
     1392            if (!isfinite(source->pixels->data.F32[iy][ix])) {
     1393                continue;
     1394            }
     1395
     1396            float fy = source->pixels->data.F32[iy][ix];
     1397            float fm = source->modelFlux->data.F32[iy][ix];
     1398            float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
     1399
     1400            YY += PS_SQR(fy) * wt;
     1401            YM += fm * fy * wt;
     1402            MM += PS_SQR(fm) * wt;
     1403        }
     1404    }
     1405
     1406    float Io = YM / MM;
     1407    PAR[PM_PAR_I0] = Io;
     1408
     1409    pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);  // NOTE : 1687 allocs in here
     1410
     1411    return true;
     1412}
     1413
     1414
Note: See TracChangeset for help on using the changeset viewer.