Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 35771) +++ trunk/psphot/src/psphotBlendFit.c (revision 36086) @@ -273,6 +273,6 @@ pmSource *source = sources->data[i]; +# if (0) int TEST_ON = false; -# if (0) # define TEST_X 653 # define TEST_Y 466 @@ -342,7 +342,9 @@ if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) { +# if (PS_TRACE_ON) if (TEST_ON) { psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); } +# endif psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); Next ++; @@ -351,7 +353,9 @@ } else { if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) { +# if (PS_TRACE_ON) if (TEST_ON) { psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); } +# endif source->type = PM_SOURCE_TYPE_STAR; psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); Index: trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- trunk/psphot/src/psphotDeblendSatstars.c (revision 35771) +++ trunk/psphot/src/psphotDeblendSatstars.c (revision 36086) @@ -128,4 +128,6 @@ // int display = psphotKapaChannel (1); // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0); + + int Nsatstar = 0; // examine sources in decreasing SN order @@ -179,4 +181,6 @@ source->mode2 |= PM_SOURCE_MODE2_SATSTAR_PROFILE; // and we have in fact subtracted the profile + Nsatstar ++; + // XXX visualize, model, and subtract // if (!psphotVisualRadialProfileSatstar (source, maskVal)) { @@ -192,5 +196,5 @@ psFree (index); - psLogMsg ("psphot", PS_LOG_INFO, "deblend satstar: %f sec\n", psTimerMark ("psphot.deblend.sat")); + psLogMsg ("psphot", PS_LOG_INFO, "deblend %d satstars: %f sec\n", Nsatstar, psTimerMark ("psphot.deblend.sat")); return true; } Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 35771) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 36086) @@ -65,4 +65,6 @@ psphotVisualShowImage(readout); + // psphotSaveImage (NULL, readout->image, "test.01.fits"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); @@ -123,7 +125,11 @@ } + // perform full extended source non-linear fits? + bool isInteractive = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_INTERACTIVE"); + if (!status) isInteractive = false; + // Define source fitting parameters for extended source fits pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); - fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->mode = PM_SOURCE_FIT_EXT_AND_SKY; fitOptions->saveCovariance = true; // XXX make this a user option? fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix @@ -134,4 +140,5 @@ fitOptions->gainFactorMode = gainFactorMode; fitOptions->chisqConvergence = chisqConvergence; + fitOptions->isInteractive = isInteractive; // maskVal is used to test for rejected pixels, and must include markVal @@ -296,41 +303,42 @@ PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail -// set this to 0 to run without threading -# if (1) - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree(AnalysisRegion); - psFree (fitOptions); - psFree (models); - psphotSersicModelClassCleanup(); - return false; - } -# else - if (!psphotExtendedSourceFits_Threaded(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree(AnalysisRegion); - psFree (fitOptions); - psFree (models); - psphotSersicModelClassCleanup(); - return false; +// XXX TEST + if (!isInteractive) { + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + psFree (fitOptions); + psFree (models); + psphotSersicModelClassCleanup(); + return false; + } + } else { + // run without threading + if (!psphotExtendedSourceFits_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + psFree (fitOptions); + psFree (models); + psphotSersicModelClassCleanup(); + return false; + } + psScalar *scalar = NULL; + scalar = job->args->data[9]; + Next += scalar->data.S32; + scalar = job->args->data[10]; + Nconvolve += scalar->data.S32; + scalar = job->args->data[11]; + NconvolvePass += scalar->data.S32; + scalar = job->args->data[12]; + Nplain += scalar->data.S32; + scalar = job->args->data[13]; + NplainPass += scalar->data.S32; + scalar = job->args->data[14]; + Nfaint += scalar->data.S32; + scalar = job->args->data[15]; + Nfail += scalar->data.S32; + psFree(job->args->data[3]); // iterator allocated above + psFree(job); } - psScalar *scalar = NULL; - scalar = job->args->data[9]; - Next += scalar->data.S32; - scalar = job->args->data[10]; - Nconvolve += scalar->data.S32; - scalar = job->args->data[11]; - NconvolvePass += scalar->data.S32; - scalar = job->args->data[12]; - Nplain += scalar->data.S32; - scalar = job->args->data[13]; - NplainPass += scalar->data.S32; - scalar = job->args->data[14]; - Nfaint += scalar->data.S32; - scalar = job->args->data[15]; - Nfail += scalar->data.S32; - psFree(job->args->data[3]); // iterator allocated above - psFree(job); -# endif } @@ -378,4 +386,6 @@ psphotSersicModelClassCleanup(); + // psphotSaveImage (NULL, readout->image, "test.02.fits"); + psphotVisualShowResidualImage (readout, false); @@ -521,4 +531,8 @@ psFree (source->modelFlux); source->modelFlux = NULL; + // if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) { + // fprintf (stderr, "test object\n"); + // } + modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal); if (!modelFit) { @@ -544,7 +558,5 @@ float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]); - if (false && (source->peak->xf > 1100) && - (source->peak->xf < 1400) && - (source->peak->yf < 245)) { + if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) { fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", source->pixels->col0, source->pixels->row0, @@ -579,5 +591,5 @@ if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue; - if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; + // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue; @@ -598,4 +610,11 @@ pmSourceCacheModel (source, maskVal); + pmModel *model = source->modelFits->data[0]; + int flags = 0xffffffff; + if (model) { + flags = model->flags; + } + + fprintf (stderr, "failed to fit extended source model to object %d @ %f, %f (%x)\n", source->id, source->moments->Mx, source->moments->My, flags); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); @@ -615,6 +634,8 @@ // cache the model flux if (source->modelEXT->isPCM) { + // fprintf (stderr, "subtract PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My); pmPCMCacheModel (source, maskVal, psfSize); } else { + // fprintf (stderr, "subtract non-PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My); pmSourceCacheModel (source, maskVal); } Index: trunk/psphot/src/psphotRadiusChecks.c =================================================================== --- trunk/psphot/src/psphotRadiusChecks.c (revision 35771) +++ trunk/psphot/src/psphotRadiusChecks.c (revision 36086) @@ -171,6 +171,9 @@ # undef MIN_WINDOW +// XXX EAM : 20130724 : for a test, double the window size parameters +// # define MIN_WINDOW 10.0 +// # define SCALE1 7.0 # define MIN_WINDOW 5.0 -# define SCALE1 5.0 +# define SCALE1 3.0 # define PAD_WINDOW 3.0 Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 35771) +++ trunk/psphot/src/psphotSourceFits.c (revision 36086) @@ -559,4 +559,10 @@ # define TIMING 0 +# define EXTRA_VERBOSE 0 + +bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize); +bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); +bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); +bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { @@ -580,5 +586,16 @@ float t1, t2, t4, t5; + t1 = t2 = t4 = t5 = 0.0; if (TIMING) { psTimerStart ("psphotFitPCM"); } + + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + options.mode = PM_SOURCE_FIT_NO_INDEX; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133 + options.mode = PM_SOURCE_FIT_EXT_AND_SKY; + // if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include + // that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image) + // thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init + } else { + options.mode = PM_SOURCE_FIT_EXT_AND_SKY; + } pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize); @@ -594,7 +611,18 @@ if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { // use the source moments, etc to guess basic model parameters - if (!psphotSersicModelClassGuessPCM (pcm, source)) { + if (!psphotSersicModelGuessPCM (pcm, source, maskVal, psfSize)) { psFree (pcm); model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } + if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); } + + // psphotFitSersicShapeAndIndex (pcm, readout, source, fitOptions, maskVal, markVal, psfSize); + options.mode = PM_SOURCE_FIT_NO_INDEX; + if (!psphotFitSersicShapeAndIndexGridAuto (pcm, readout, source, &options, maskVal, markVal, psfSize)) { + psFree (pcm); + model->flags |= PM_MODEL_STATUS_BADARGS; + psError(PS_ERR_UNKNOWN, true, "Failed to find a index & shape"); + psErrorClear (); // clear the polynomial error return model; } @@ -608,13 +636,4 @@ } - if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); } - - if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { - options.mode = PM_SOURCE_FIT_NO_INDEX; - } else { - options.mode = PM_SOURCE_FIT_EXT; - } - // update the pcm elements if we have changed the circumstance (options.mode or source->pixels) - pmPCMupdate(pcm, source, &options, model); if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); } @@ -629,6 +648,10 @@ 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); } - - // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); + if (EXTRA_VERBOSE && !TIMING) { + int nPixBig = source->pixels->numCols * source->pixels->numRows; + float *PAR = model->params->data.F32; + fprintf (stderr, "source %d : %f - %f %f - %f %f %f - %f | nIter: %2d, radius: %6.1f, npix: %5d of %5d, chisq %f\n", source->id, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1], model->nIter, model->fitRadius, model->nPix, nPixBig, model->chisqNorm); + } + psFree (pcm); @@ -748,2 +771,472 @@ return true; } + +// float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083}; +float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125}; +# define N_INDEX_GUESS_INV 7 + +// float reffGuess[] = {3.0, 10.0, 20.0, 30.0, 40.0}; +float reffGuess[] = {0.5, 0.75, 1.0, 1.4, 2.0}; +# define N_REFF_GUESS 5 + +// A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index +// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations +bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) { + + // we get a reasonable guess from: + // * Reff = Kron R1 + // * Rmajor / Rminor & Theta from moments + // * Io from total Kron flux + + // the guesses are used to fill in PAR: + psF32 *PAR = pcm->modelConv->params->data.F32; + + // convert the moments to Major,Minor,Theta + psEllipseMoments moments; + + if (!isfinite(source->moments->Mrf)) return false; + if (!isfinite(source->moments->Mxx)) return false; + if (!isfinite(source->moments->Mxy)) return false; + if (!isfinite(source->moments->Myy)) return false; + + moments.x2 = source->moments->Mxx; + moments.y2 = source->moments->Myy; + moments.xy = source->moments->Mxy; + + // limit axis ratio < 20.0 + psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0); + + if (0) { + psEllipseAxes guessAxes; + guessAxes.major = source->moments->Mrf; + guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major; + guessAxes.theta = momentAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + // sky is zero (no longer fitted, but not yet deprecated) + PAR[PM_PAR_SKY] = 0.0; + + // for the index loop, use Io = 1.0, use fitted values to determine Io + PAR[PM_PAR_I0] = 1.0; + + float xMin = NAN; + float iMin = NAN; + float sMin = NAN; + float rMin = NAN; + + // loop over index and Reff, keeping the ARatio and Theta constant? + // loop over index guesses and find the best fit + for (int j = 0; j < N_REFF_GUESS; j++) { + + psEllipseAxes guessAxes; + guessAxes.major = reffGuess[j] * source->moments->Mrf; + guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major); + guessAxes.theta = momentAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + + for (int i = 0; i < N_INDEX_GUESS_INV; i++) { + PAR[PM_PAR_7] = indexGuessInv[i]; + + // generated the modelFlux + pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize); + + float YY = 0.0; + float YM = 0.0; + float MM = 0.0; + bool usePoisson = false; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[iy][ix] == 0) { + continue; + } + // skip nan value points + if (!isfinite(source->pixels->data.F32[iy][ix])) { + continue; + } + + float fy = source->pixels->data.F32[iy][ix]; + float fm = source->modelFlux->data.F32[iy][ix]; + float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0; + + YY += PS_SQR(fy) * wt; + YM += fm * fy * wt; + MM += PS_SQR(fm) * wt; + } + } + + float Io = YM / MM; + float Chisq = YY - 2 * Io * YM + Io * Io * MM; + if (isnan(xMin) || (Chisq < xMin)) { + xMin = Chisq; + iMin = Io; + sMin = indexGuessInv[i]; + rMin = reffGuess[j]; + } + // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin); + // fprintf (stderr, "\n"); + } + } + + { + psEllipseAxes guessAxes; + guessAxes.major = rMin * source->moments->Mrf; + guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major); + guessAxes.theta = momentAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + } + + PAR[PM_PAR_I0] = iMin; + PAR[PM_PAR_7] = sMin; + + return true; +} + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + for (int i = 0; i < 3; i++) { + // fit EXT (not PSF) model (set/unset the pixel mask) + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 2; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) { + float *PAR = model->params->data.F32; + 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]); + } + + // fit EXT (not PSF) model (set/unset the pixel mask) + options.mode = PM_SOURCE_FIT_INDEX; + // options.mode = PM_SOURCE_FIT_EXT_AND_SKY; + options.nIter = 30; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) { + float *PAR = model->params->data.F32; + 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]); + } + } + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, fitOptions, model); + + return true; +} + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + psF32 *PAR = pcm->modelConv->params->data.F32; + + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 7; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + // we have been provided a guess at the index (P[7]) from the list of indexGuessInv + + // find the matching indexGuessInv + int nStart = -1; + for (int i = 0; i < N_INDEX_GUESS_INV; i++) { + if (fabs(PAR[PM_PAR_7] - indexGuessInv[i]) < 0.01) { + nStart = i; + break; + } + } + if (nStart == -1) { + fprintf (stderr, "WARNING: could not find start guess %f\n", PAR[PM_PAR_7]); + return false; + } + + psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32); + psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32); + + PAR[PM_PAR_7] = indexGuessInv[nStart]; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.1; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + if (chi2->data.F32[1] < chi2->data.F32[2]) { + if (nStart == N_INDEX_GUESS_INV - 1) { + PAR[PM_PAR_7] = 0.11; + } else { + PAR[PM_PAR_7] = indexGuessInv[nStart + 1]; + } + } else { + if (nStart == 0) { + PAR[PM_PAR_7] = 0.52; + } else { + PAR[PM_PAR_7] = indexGuessInv[nStart - 1]; + } + } + + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); + if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola"); + psFree (chi2); + psFree (Sidx); + psFree (poly); + return false; + } + + // where is the minimum of this polynomial fit? + float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0; + + // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev) + Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1); + PAR[PM_PAR_7] = Smin; + + // return to the original fitting mode (fitOptions) + pmPCMupdate(pcm, source, fitOptions, model); + + psFree (chi2); + psFree (Sidx); + psFree (poly); + + return true; +} + + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndexGridAutoScaled (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + psF32 *PAR = pcm->modelConv->params->data.F32; + + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 5; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + float parStart[8]; + for (int i = 0; i < 8; i++) parStart[i] = PAR[i]; + + // we start with a guess at the index (P[7]) + + // get chisq for P[7], P[7]*1.1, P[7]*1.25 (or *0.75 depending on the result of 1.1) + + psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32); + psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32); + + PAR[PM_PAR_7] = parStart[7]; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + float fI = 1.1; + PAR[PM_PAR_7] = parStart[7]*fI; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + if (chi2->data.F32[1] < chi2->data.F32[0]) { + fI = 1.3; + } else { + fI = 1.0 / 1.3; + } + + PAR[PM_PAR_7] = parStart[7]*fI; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + // can we fit the 3 pts with a parabola? + int nTry = 0; + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); + while (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) { + psErrorClear (); // clear the polynomial error + if (nTry > 4) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola"); + psFree (chi2); + psFree (Sidx); + psFree (poly); + return false; + } + fI = (fI < 1.0) ? fI / 1.3 : fI * 1.3; + PAR[PM_PAR_7] = parStart[7]*fI; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + nTry ++; + } + + // where is the minimum of this polynomial fit? + float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0; + + // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev) + Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1); + PAR[PM_PAR_7] = Smin; + + // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + // 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]); + + //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i]; + //// + //// for (float fI = 0.0; fI < 0.15; fI += 0.01) { + //// PAR[PM_PAR_7] = parStart[7] - fI; + //// + //// // fit EXT (not PSF) model (set/unset the pixel mask) + //// + //// pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + //// if (TIMING) { + //// float *PAR = model->params->data.F32; + //// fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]); + //// } + //// } + + // return to the original fitting mode (fitOptions) + pmPCMupdate(pcm, source, fitOptions, model); + + psFree (chi2); + psFree (Sidx); + psFree (poly); + + return true; +} + + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + psF32 *PAR = pcm->modelConv->params->data.F32; + + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 10; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32); + psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32); + + float par7[] = {0.100, 0.125, 0.150, 0.175, 0.200, 0.225, 0.250}; + for (int i = 0; i < 7; i++) { + PAR[PM_PAR_7] = par7[i]; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + psVectorAppend (Sidx, PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + } + + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); + if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola"); + psFree (chi2); + psFree (Sidx); + psFree (poly); + return false; + } + + // where is the minimum of this polynomial fit? + fprintf (stderr, "fit1d: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + float Smin = -0.5 * poly->coeff[1] / poly->coeff[2]; + + // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev) + Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1); + PAR[PM_PAR_7] = Smin; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + 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]); + + //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i]; + //// + //// for (float fI = 0.0; fI < 0.15; fI += 0.01) { + //// PAR[PM_PAR_7] = parStart[7] - fI; + //// + //// // fit EXT (not PSF) model (set/unset the pixel mask) + //// + //// pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + //// if (TIMING) { + //// float *PAR = model->params->data.F32; + //// fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]); + //// } + //// } + + // return to the original fitting mode (fitOptions) + pmPCMupdate(pcm, source, fitOptions, model); + + psFree (chi2); + psFree (Sidx); + psFree (poly); + + return true; +} + +