Index: /trunk/psphot/src/models/pmModel_STRAIL.c =================================================================== --- /trunk/psphot/src/models/pmModel_STRAIL.c (revision 35558) +++ /trunk/psphot/src/models/pmModel_STRAIL.c (revision 35559) @@ -475,5 +475,5 @@ //fixed I think...no good way of guessing as far as I can tell -bool PM_MODEL_GUESS (pmModel *model, pmSource *source) +bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { pmMoments *Smoments = source->moments; Index: /trunk/psphot/src/models/pmModel_TEST1.c =================================================================== --- /trunk/psphot/src/models/pmModel_TEST1.c (revision 35558) +++ /trunk/psphot/src/models/pmModel_TEST1.c (revision 35559) @@ -123,5 +123,5 @@ // make an initial guess for parameters -bool PM_MODEL_GUESS (pmModel *model, pmSource *source) +bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { pmMoments *moments = source->moments; Index: /trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceFits.c (revision 35558) +++ /trunk/psphot/src/psphotExtendedSourceFits.c (revision 35559) @@ -90,4 +90,21 @@ assert (markVal); + // source fitting parameters for extended source fits + int fitIter = psMetadataLookupS32(&status, recipe, "EXT_FIT_ITER"); // Max number of fit iterations + assert (status && fitIter > 0); + + float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance + if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) { + fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance + if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) { + psAbort("PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive"); + } + } + + float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance + if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) { + fitMaxTol = 1.0; + } + // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; @@ -181,4 +198,8 @@ PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, fitIter, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, fitMinTol, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, fitMaxTol, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); @@ -293,6 +314,9 @@ psRegion *region = job->args->data[4]; int psfSize = PS_SCALAR_VALUE(job->args->data[5],S32); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); + int fitIter = PS_SCALAR_VALUE(job->args->data[6],S32); + float fitMinTol = PS_SCALAR_VALUE(job->args->data[7],F32); + float fitMaxTol = PS_SCALAR_VALUE(job->args->data[8],F32); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[9],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[10],PS_TYPE_IMAGE_MASK_DATA); // Define source fitting parameters for extended source fits @@ -301,10 +325,7 @@ fitOptions->saveCovariance = true; // XXX make this a user option? fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix - - // XXX for now, use the defaults for the rest: - // fitOptions->nIter = fitIter; - // fitOptions->tol = fitTol; - // fitOptions->poissonErrors = poisson; - // fitOptions->weight = PS_SQR(skySig); + fitOptions->nIter = fitIter; + fitOptions->minTol = fitMinTol; + fitOptions->maxTol = fitMaxTol; // choose the sources of interest @@ -335,4 +356,7 @@ // set the fit radius based on the first radial moment (also sets the mask pixels) psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs + + // XXX WATCH OUT HERE!! + // fitRadius = 30; // UPDATE : we have changed the moments calculation. There is now an iteration within @@ -498,23 +522,23 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[8]; + scalar = job->args->data[11]; scalar->data.S32 = Next; - scalar = job->args->data[9]; + scalar = job->args->data[12]; scalar->data.S32 = Nconvolve; - scalar = job->args->data[10]; + scalar = job->args->data[13]; scalar->data.S32 = NconvolvePass; - scalar = job->args->data[11]; + scalar = job->args->data[14]; scalar->data.S32 = Nplain; - scalar = job->args->data[12]; + scalar = job->args->data[15]; scalar->data.S32 = NplainPass; - scalar = job->args->data[13]; + scalar = job->args->data[16]; scalar->data.S32 = Nfaint; - scalar = job->args->data[14]; + scalar = job->args->data[17]; scalar->data.S32 = Nfail; Index: /trunk/psphot/src/psphotFakeSources.c =================================================================== --- /trunk/psphot/src/psphotFakeSources.c (revision 35558) +++ /trunk/psphot/src/psphotFakeSources.c (revision 35559) @@ -24,5 +24,5 @@ pmModelType modelType = pmModelClassGetType ("PS_MODEL_QGAUSS"); - source->modelPSF = pmSourceModelGuess (source, modelType); + source->modelPSF = pmSourceModelGuess (source, modelType, 0, 0); sources->data[i] = source; } Index: /trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- /trunk/psphot/src/psphotMakeResiduals.c (revision 35558) +++ /trunk/psphot/src/psphotMakeResiduals.c (revision 35559) @@ -105,7 +105,7 @@ if (model == NULL) continue; // model must be defined - psImage *image = psImageCopy (NULL, source->pixels, PS_TYPE_F32); - psImage *variance = psImageCopy (NULL, source->variance, PS_TYPE_F32); - psImage *mask = psImageCopy (NULL, source->maskView, PS_TYPE_IMAGE_MASK); + psImage *image = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + psImage *mask = psImageCopy (NULL, source->maskView ? source->maskView : source->maskObj, PS_TYPE_IMAGE_MASK); + psImage *variance = psImageCopy (NULL, source->variance ? source->variance : source->pixels, PS_TYPE_F32); pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal); Index: /trunk/psphot/src/psphotModelTestArguments.c =================================================================== --- /trunk/psphot/src/psphotModelTestArguments.c (revision 35558) +++ /trunk/psphot/src/psphotModelTestArguments.c (revision 35559) @@ -151,4 +151,10 @@ } + // chip selection is used to limit chips to be processed + if ((N = psArgumentGet (argc, argv, "-moments-radius"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "TEST_MOMENTS_RADIUS", PS_DATA_BOOL, "", true); + } + // if these command-line options are supplied, load the file name lists into config->arguments // override any configuration-specified source for these files @@ -160,4 +166,6 @@ pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC", "-src", "-srclist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "SRCTEXT", "-srctext", "-srctextlist"); + + pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); if (argc == 1) { Index: /trunk/psphot/src/psphotModelTestReadout.c =================================================================== --- /trunk/psphot/src/psphotModelTestReadout.c (revision 35558) +++ /trunk/psphot/src/psphotModelTestReadout.c (revision 35559) @@ -5,4 +5,6 @@ bool status; + + psTimerStart ("modelTest"); pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); @@ -41,4 +43,11 @@ } + // load the psf model, if suppled. FWHM_MAJ,FWHM_MIN,etc are determined and saved on + // readout->analysis. NOTE: this function currently only loads from PSPHOT.PSF.LOAD + if (!psphotLoadPSF (config, view, filerule)) { // ??? need to supply 2 ? + psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); + return psphotReadoutCleanup (config, view, filerule); + } + float MIN_KRON_RADIUS = 5.0; @@ -53,4 +62,8 @@ // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); // find the currently selected readout @@ -84,5 +97,5 @@ // define the source of interest - float xObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X"); + float xObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X"); float yObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y"); if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined"); @@ -99,9 +112,12 @@ source->type = PM_SOURCE_TYPE_EXTENDED; + bool TEST_FIT_CONVOLVED = psMetadataLookupBool (&status, recipe, "TEST_FIT_CONVOLVED"); + // find the local sky status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); if (!status) psAbort("pmSourceLocalSky error"); - { + bool testMomentsRadius = psMetadataLookupBool (&status, config->arguments, "TEST_MOMENTS_RADIUS"); + if (testMomentsRadius) { // XXX I want to test an iterative aperture for brighter sources float radius = mRADIUS; @@ -109,5 +125,5 @@ // get the source moments - status = pmSourceMoments (source, radius, 0.0, 0.0, MIN_KRON_RADIUS, maskVal); + status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal); if (!status) psAbort("psSourceMoments error"); @@ -121,5 +137,5 @@ // get the source moments - status = pmSourceMoments (source, mRADIUS, 0.0, 0.0, MIN_KRON_RADIUS, maskVal); + status = pmSourceMoments (source, mRADIUS, 0.25*mRADIUS, 0.0, MIN_KRON_RADIUS, maskVal); if (!status) psAbort("psSourceMoments error"); @@ -135,6 +151,14 @@ fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + if (psf) { + // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo) + source->peak->rawFlux = source->moments->Peak; + float Io = source->moments->Peak; + source->modelPSF = pmModelFromPSFforXY(psf, xObj, yObj, Io); + } + // get the initial model parameter guess - pmModel *model = pmSourceModelGuess (source, modelType); + pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal); source->modelEXT = model; @@ -159,12 +183,15 @@ fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y); - // list model input shape - psEllipseShape shape; - shape.sx = 1.4 / model->params->data.F32[4]; - shape.sy = 1.4 / model->params->data.F32[5]; - shape.sxy = model->params->data.F32[6]; - axes = psEllipseShapeToAxes (shape, 20.0); - - fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) { + fprintf (stderr, "guess: %f @ (%f, %f)\n", params[6]*180/M_PI, params[4], params[5]); + } else { + // list model input shape + psEllipseShape shape; + shape.sx = 1.4 / model->params->data.F32[4]; + shape.sy = 1.4 / model->params->data.F32[5]; + shape.sxy = model->params->data.F32[6]; + axes = psEllipseShapeToAxes (shape, 20.0); + fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + } fprintf (stderr, "input parameters: \n"); @@ -190,5 +217,24 @@ fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix - status = pmSourceFitModel (source, model, fitOptions, maskVal); + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + fitOptions->mode = PM_SOURCE_FIT_NO_INDEX; + } + if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) { + fitOptions->mode = PM_SOURCE_FIT_TRAIL; + } + + if (TEST_FIT_CONVOLVED) { + pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + psphotSersicModelClassGuessPCM (pcm, source); + } else { + pmSourceModelGuessPCM (pcm, source, maskVal, markVal); + } + pmPCMupdate(pcm, source, fitOptions, model); + pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize); + psFree (pcm); + } else { + status = pmSourceFitModel (source, model, fitOptions, maskVal); + } // measure the source mags Index: /trunk/psphot/src/psphotSetThreads.c =================================================================== --- /trunk/psphot/src/psphotSetThreads.c (revision 35558) +++ /trunk/psphot/src/psphotSetThreads.c (revision 35559) @@ -40,5 +40,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 15); + task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 18); task->function = &psphotExtendedSourceFits_Threaded; psThreadTaskAdd(task); Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 35558) +++ /trunk/psphot/src/psphotSourceFits.c (revision 35559) @@ -457,5 +457,5 @@ // use the source moments, etc to guess basic model parameters - pmModel *model = pmSourceModelGuess (source, modelType); + pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal); if (!model) { psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); @@ -625,5 +625,5 @@ model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i]; - if (!model->modelGuess(model, source)) { + if (!model->modelGuess(model, source, maskVal, markVal)) { model->flags |= PM_MODEL_STATUS_BADARGS; return false; @@ -648,5 +648,5 @@ model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin]; - model->modelGuess(model, source); + model->modelGuess(model, source, maskVal, markVal); return true; @@ -677,5 +677,5 @@ model->params->data.F32[PM_PAR_7] = indexGuess[i]; - if (!model->modelGuess(model, source)) { + if (!model->modelGuess(model, source, maskVal, markVal)) { model->flags |= PM_MODEL_STATUS_BADARGS; return false;