Index: /branches/eam_branches/ipp-20130904/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36217) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36218) @@ -421,4 +421,13 @@ } psphotStackOptions; +typedef struct { + float fRmajorMin; + float fRmajorMax; + float fRmajorDel; + float fRminorMin; + float fRminorMax; + float fRminorDel; +} psphotGalaxyShapeOptions; + /*** psphotStackMatchPSF prototypes ***/ bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view); @@ -517,6 +526,7 @@ bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf); bool psphotGalaxyShape_Threaded (psThreadJob *job); -bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize); +bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize); bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize); +psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc(); bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule); Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c (revision 36217) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c (revision 36218) @@ -83,5 +83,15 @@ maskVal |= markVal; + // bit-mask to mark pixels not used in analysis + psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc(); + opt->fRmajorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MIN"); if (!status) opt->fRmajorMin = 0.5; + opt->fRmajorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MAX"); if (!status) opt->fRmajorMax = 2.0; + opt->fRmajorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_DEL"); if (!status) opt->fRmajorDel = 0.1; + opt->fRminorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MIN"); if (!status) opt->fRminorMin = 0.5; + opt->fRminorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MAX"); if (!status) opt->fRminorMax = 2.0; + opt->fRminorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_DEL"); if (!status) opt->fRminorDel = 0.1; + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + // NOTE: this is only used if we are NOT smoothing with a 1D Gaussian int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); assert (status); @@ -119,4 +129,5 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources psArrayAdd(job->args, 1, fitOptions); + psArrayAdd(job->args, 1, opt); PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); @@ -156,4 +167,5 @@ psFree (cellGroups); psFree (fitOptions); + psFree (opt); psLogMsg ("psphot.galaxy", PS_LOG_WARN, "measure galaxy shapes : %f sec for %ld objects\n", psTimerMark ("psphot.galaxy"), sources->n); @@ -166,7 +178,8 @@ psArray *sources = job->args->data[1]; pmSourceFitOptions *fitOptions = job->args->data[2]; - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); - int psfSize = PS_SCALAR_VALUE(job->args->data[5],S32); + psphotGalaxyShapeOptions *opt = job->args->data[3]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + int psfSize = PS_SCALAR_VALUE(job->args->data[6],S32); float fitRadius; @@ -202,5 +215,6 @@ // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr // do the following for a set of shapes (Ex,Ey) - psphotGalaxyShapeGrid (source, fitOptions, maskVal, psfSize); + psphotGalaxyShapeGrid (source, fitOptions, opt, maskVal, psfSize); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); @@ -211,5 +225,5 @@ } -bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize) { +bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) { // we use the 0th model to define the initial guess shape @@ -230,28 +244,23 @@ // I have some source guess (e0, e1, e2) psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType); - psEllipsePol guessPol = psEllipseAxesToPol (guessAxes); - - // I am going to retain the value of e0, and grid search around e1,e2 - // I need to study a valid range (and distribution?) for e1,e2) if (!source->galaxyFits) { - source->galaxyFits = pmSourceGalaxyFitsAlloc(); - } - - for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) { - for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) { + source->galaxyFits = pmSourceGalaxyFitsAlloc(); + } + + for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) { + for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) { - psEllipsePol testPol = guessPol; - testPol.e1 += de1; - testPol.e2 += de2; - - psEllipseAxes axes = psEllipsePolToAxes (testPol, 0.1); - - pmPSF_AxesToModel (PAR, axes, modelType); - + psEllipseAxes testAxes = guessAxes; + testAxes.major = guessAxes.major * fRmajor; + testAxes.minor = guessAxes.minor * fRminor; + + pmPSF_AxesToModel (PAR, testAxes, modelType); + // where do the results go?? psphotGalaxyShapeSource (pcm, source, maskVal, psfSize); } } + psFree (pcm); return true; @@ -267,5 +276,7 @@ // generate the modelFlux + // need to reset here each time since we assign below to calculate the flux pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize); + pcm->modelConv->params->data.F32[PM_PAR_I0] = 1.0; int nPix = 0; @@ -273,5 +284,5 @@ float YM = 0.0; float MM = 0.0; - bool usePoisson = false; + bool usePoisson = true; for (int iy = 0; iy < source->pixels->numRows; iy++) { @@ -306,8 +317,10 @@ float Chisq = (YY - 2 * Io * YM + Io * Io * MM) / (float) nPix; - // fprintf (stderr, "Io: %f, Chisq: %f\n", Io, Chisq); - - psVectorAppend (source->galaxyFits->Io, Io); - psVectorAppend (source->galaxyFits->dIo, dIo); + pcm->modelConv->params->data.F32[PM_PAR_I0] = Io; + float flux = pcm->modelConv->modelFlux (pcm->modelConv->params); + float dflux = flux * (dIo / Io); + + psVectorAppend (source->galaxyFits->Flux, flux); + psVectorAppend (source->galaxyFits->dFlux, dflux); psVectorAppend (source->galaxyFits->chisq, Chisq); source->galaxyFits->nPix = nPix; @@ -315,2 +328,24 @@ return true; } + +/**** support functions ****/ + +void psphotGalaxyShapeOptionsFree (psphotGalaxyShapeOptions *opt) { + return; +} + +psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc() +{ + psphotGalaxyShapeOptions *opt = (psphotGalaxyShapeOptions *) psAlloc(sizeof(psphotGalaxyShapeOptions)); + psMemSetDeallocator(opt, (psFreeFunc) psphotGalaxyShapeOptionsFree); + + opt->fRmajorMin = 0.5; + opt->fRmajorMax = 2.0; + opt->fRmajorDel = 0.1; + opt->fRminorMin = 0.5; + opt->fRminorMax = 2.0; + opt->fRminorDel = 0.1; + + return opt; +} + Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotSetThreads.c (revision 36217) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotSetThreads.c (revision 36218) @@ -57,5 +57,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPES", 6); + task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPES", 7); task->function = &psphotGalaxyShape_Threaded; psThreadTaskAdd(task);