IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36318


Ignore:
Timestamp:
Nov 28, 2013, 4:09:45 AM (13 years ago)
Author:
eugene
Message:

make sure the convolution scale is passed to all functions that need it; update sersic search a bit

Location:
branches/eam_branches/ipp-20130904/psphot/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c

    r36315 r36318  
    120120    }
    121121
     122    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio
     123    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
     124        fitNsigmaConv = 5.0;
     125    }
     126
    122127    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
    123128    if (!status) {
    … …  
    144149    fitOptions->minTol         = fitMinTol;
    145150    fitOptions->maxTol         = fitMaxTol;
     151    fitOptions->nsigma         = fitNsigmaConv;
    146152
    147153    fitOptions->gainFactorMode   = gainFactorMode;
    … …  
    656662        if (source->modelEXT->isPCM) {
    657663            // fprintf (stderr, "subtract PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);
    658             pmPCMCacheModel (source, maskVal, psfSize);
     664            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
    659665        } else {
    660666            // fprintf (stderr, "subtract non-PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);

  • branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c

    r36315 r36318  
    9797    assert (status);
    9898
     99    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio
     100    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
     101        fitNsigmaConv = 5.0;
     102    }
     103
    99104    // Define source fitting parameters for extended source fits
    100105    // we are not doing LMM fitting, so most options are irrelevant
    … …  
    102107    fitOptions->mode           = PM_SOURCE_FIT_EXT_AND_SKY;
    103108    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
     109    fitOptions->nsigma         = fitNsigmaConv;
    104110
    105111    // Poisson or Constant weight for chisq tests?
    … …  
    298304
    299305        // cache the model flux
    300         pmPCMCacheModel (source, maskVal, psfSize);
     306        pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
    301307    }
    302308
    … …  
    316322    // need to reset here each time since we assign below to calculate the flux
    317323    // XXX note that this does not add sky to model
    318     pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
     324    pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
    319325    pcm->modelConv->params->data.F32[PM_PAR_I0] = 1.0;
    320326       

  • branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c

    r36315 r36318  
    653653    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
    654654
    655     psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
    656     if (pcm->modelConv->nIter == fitOptions->nIter) {
    657         psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
    658     }
     655    // psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
     656    // if (pcm->modelConv->nIter == fitOptions->nIter) {
     657    // psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
     658    // }
    659659    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
    660660
    … …  
    787787}
    788788
    789 // float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
    790 float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
    791 # define N_INDEX_GUESS_INV 7
    792 
    793 // float reffGuess[] = {3.0, 10.0, 20.0, 30.0, 40.0};
    794 float reffGuess[] = {0.5, 0.75, 1.0, 1.4, 2.0};
     789// 0.5 / n for (1.0, 1.25, 1.66, 2.0, 3.33, 4.0)
     790// float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
     791
     792// 0.5 / n for (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0)
     793float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 1.666, 0.125, 0.10, 0.0833};
     794float indexGuessR1q[] = {1.06, 1.19, 1.335, 1.48, 1.840, 2.290, 2.84, 3.5300};
     795# define N_INDEX_GUESS_INV 8
     796
     797// we are going to guess in fractions about the R1-based guess
     798float reffGuess[] = {0.8, 0.9, 1.0, 1.12, 1.25};
    795799# define N_REFF_GUESS 5
    796800
    … …  
    800804
    801805    // we get a reasonable guess from:
    802     // * Reff = Kron R1
     806    // * Reff = Kron R1 / Q(index) -- Q comes from Graham & Driver
    803807    // * Rmajor / Rminor & Theta from moments
    804808    // * Io from total Kron flux
    … …  
    822826    psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0);
    823827
    824     if (0) {
    825         psEllipseAxes guessAxes;
    826         guessAxes.major = source->moments->Mrf;
    827         guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
    828         guessAxes.theta = momentAxes.theta;
    829 
    830         if (!isfinite(guessAxes.major)) return false;
    831         if (!isfinite(guessAxes.minor)) return false;
    832         if (!isfinite(guessAxes.theta)) return false;
    833 
    834         // convert the major,minor,theta to shape parameters for an Reff-like model
    835         pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
    836     }
    837 
    838828    // set the model position
    839829    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
    … …  
    855845    // loop over index guesses and find the best fit
    856846    for (int j = 0; j < N_REFF_GUESS; j++) {
    857 
    858         psEllipseAxes guessAxes;
    859         guessAxes.major = reffGuess[j] * source->moments->Mrf;
    860         guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
    861         guessAxes.theta = momentAxes.theta;
    862 
    863         if (!isfinite(guessAxes.major)) return false;
    864         if (!isfinite(guessAxes.minor)) return false;
    865         if (!isfinite(guessAxes.theta)) return false;
    866 
    867         // convert the major,minor,theta to shape parameters for an Reff-like model
    868         pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
    869 
    870847        for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
    871848            PAR[PM_PAR_7] = indexGuessInv[i];
    872849
     850            psEllipseAxes guessAxes;
     851            guessAxes.major = reffGuess[j] * source->moments->Mrf / indexGuessR1q[i];
     852            guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
     853            guessAxes.theta = momentAxes.theta;
     854
     855            if (!isfinite(guessAxes.major)) return false;
     856            if (!isfinite(guessAxes.minor)) return false;
     857            if (!isfinite(guessAxes.theta)) return false;
     858
     859            // convert the major,minor,theta to shape parameters for an Reff-like model
     860            pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
     861
    873862            // generated the modelFlux
    874863            // XXX note that this does not add sky to model
    875             pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
     864            pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
    876865       
    877866            float YY = 0.0;
    … …  
    911900                iMin = Io;
    912901                sMin = indexGuessInv[i];
    913                 rMin = reffGuess[j];
     902                rMin = reffGuess[j] / indexGuessR1q[i];
    914903            }
    915904            // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
    … …  
    10231012    psVectorAppend (chi2, model->chisqNorm);
    10241013
    1025     PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.1;
     1014    PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : indexGuessInv[nStart - 1];
    10261015    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    10271016    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]);
    … …  
    10291018    psVectorAppend (chi2, model->chisqNorm);
    10301019
    1031     PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55;
     1020    PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : indexGuessInv[nStart + 1];
    10321021    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    10331022    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]);
    10341023    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
    10351024    psVectorAppend (chi2, model->chisqNorm);
    1036 
    1037     if (chi2->data.F32[1] < chi2->data.F32[2]) {
    1038       if (nStart == N_INDEX_GUESS_INV - 1) {
    1039         PAR[PM_PAR_7] = 0.11;
    1040       } else {
    1041         PAR[PM_PAR_7] = indexGuessInv[nStart + 1];
    1042       }
    1043     } else {
    1044       if (nStart == 0) {
    1045         PAR[PM_PAR_7] = 0.52;
    1046       } else {
    1047         PAR[PM_PAR_7] = indexGuessInv[nStart - 1];
    1048       }
    1049     }
    10501025
    10511026    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
    … …  
    12921267        // generated the modelFlux
    12931268            // XXX note that this does not add sky to model
    1294         pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
     1269        pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
    12951270       
    12961271        float YY = 0.0;
    … …  
    13311306            rMin = dref;
    13321307        }
    1333         // fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
     1308        fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
    13341309    }
    13351310
    … …  
    13751350    // generated the modelFlux
    13761351            // XXX note that this does not add sky to model
    1377     pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
     1352    pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
    13781353       
    13791354    float YY = 0.0;
Note: See TracChangeset for help on using the changeset viewer.