Index: branches/eam_branches/ipp-20130904/psphot/src/Makefile.am =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/Makefile.am (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/Makefile.am (revision 36148) @@ -25,5 +25,5 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotForced psphotMinimal psphotMakePSF psphotStack psphotModelTest +bin_PROGRAMS = psphot psphotForced psphotFullForce psphotMinimal psphotMakePSF psphotStack psphotModelTest # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy @@ -35,4 +35,8 @@ psphotForced_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotForced_LDADD = libpsphot.la + +psphotFullForce_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotFullForce_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotFullForce_LDADD = libpsphot.la psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -77,4 +81,13 @@ psphotForced.c \ psphotForcedArguments.c \ + psphotParseCamera.c \ + psphotImageLoop.c \ + psphotMosaicChip.c \ + psphotCleanup.c + +# generalized forced photometry (including Kron, Petro, and Models) of specified positions given a specified psf +psphotFullForce_SOURCES = \ + psphotFullForce.c \ + psphotFullForceArguments.c \ psphotParseCamera.c \ psphotImageLoop.c \ @@ -155,4 +168,5 @@ psphotReadoutMinimal.c \ psphotForcedReadout.c \ + psphotFullForceReadout.c \ psphotMakePSFReadout.c \ psphotModelBackground.c \ @@ -198,4 +212,7 @@ psphotRadialPlot.c \ psphotKronIterate.c \ + psphotKronFlux.c \ + psphotPetroFlux.c \ + psphotGalaxyShape.c \ psphotRadialProfileWings.c \ psphotDeblendSatstars.c \ Index: branches/eam_branches/ipp-20130904/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36148) @@ -342,4 +342,7 @@ bool psphotForcedReadout(pmConfig *config, const pmFPAview *view, const char *filerule); +pmConfig *psphotFullForceArguments(int argc, char **argv); +bool psphotFullForceReadout(pmConfig *config, const pmFPAview *view, const char *filerule); + pmConfig *psphotMinimalArguments(int argc, char **argv); bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view, const char *filerule); @@ -498,4 +501,20 @@ bool psphotKronIterate_Threaded (psThreadJob *job); +bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources); +bool psphotKronFlux_Threaded (psThreadJob *job); +bool psphotKronFluxSource (pmSource *source, psImageMaskType maskVal); + +bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources); +bool psphotPetroFlux_Threaded (psThreadJob *job); +bool psphotPetroFluxSource (pmSource *source, psImageMaskType maskVal); + +bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule); +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 psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize); + bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources); Index: branches/eam_branches/ipp-20130904/psphot/src/psphotBlendFit.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotBlendFit.c (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotBlendFit.c (revision 36148) @@ -273,9 +273,10 @@ pmSource *source = sources->data[i]; -# if (0) - int TEST_ON = false; +# if (PS_TRACE_ON) # define TEST_X 653 # define TEST_Y 466 - if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { +# define TESTING 1 + int TEST_ON = false; + if (TESTING && (fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { fprintf (stderr, "test object\n"); psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); Index: branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceReadout.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceReadout.c (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceReadout.c (revision 36148) @@ -57,5 +57,9 @@ // generate a psf model for any readouts which need one - psphotFullForcePSF (config, view, filerule); + // psphotFullForcePSF (config, view, filerule); + // XXX I'm not sure we need a different algorithm here or not : we are supplying the + // sources; we should mark with a flag bit the ones we actually want to use as PSF + // stars (this means we need to supply this info in the load). + psphotChoosePSF (config, view, filerule, true); // Construct an initial model for each object, set the radius to fitRadius, set circular @@ -73,5 +77,5 @@ // measure kron fluxes - psphotKronFlux (config, view, filerule, false); + psphotKronFlux (config, view, filerule); // measure petro fluxes Index: branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c (revision 36148) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) +bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -38,5 +38,5 @@ // psAssert (psf, "missing psf?"); - if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) { + if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; @@ -46,5 +46,5 @@ } -bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) { +bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf) { bool status = false; @@ -55,5 +55,5 @@ } - psTimerStart ("psphot.kron"); + psTimerStart ("psphot.galaxy"); // determine the number of allowed threads @@ -63,14 +63,4 @@ } - bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT"); - if (!status) { - KRON_APPLY_WEIGHT = true; - } - - bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW"); - if (!status) { - KRON_APPLY_WINDOW = false; - } - // bit-masks to test for good/bad pixels psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); @@ -83,4 +73,54 @@ // 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); + + // 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; + } + + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + chisqConvergence = true; + } + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + gainFactorMode = 0; + } + + // 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_AND_SKY; + fitOptions->saveCovariance = true; // XXX make this a user option? + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + fitOptions->nIter = fitIter; + fitOptions->minTol = fitMinTol; + fitOptions->maxTol = fitMaxTol; + + fitOptions->gainFactorMode = gainFactorMode; + fitOptions->chisqConvergence = chisqConvergence; + fitOptions->isInteractive = isInteractive; // source analysis is done in S/N order (brightest first) @@ -105,11 +145,12 @@ // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE"); + psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_SHAPES"); psArrayAdd(job->args, 1, readout); psArrayAdd(job->args, 1, cells->data[j]); // sources - PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, pass, PS_TYPE_S32); + psArrayAdd(job->args, 1, fitOptions); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32); // set this to 0 to run without threading @@ -146,5 +187,5 @@ psFree (cellGroups); - psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); + psLogMsg ("psphot.galaxy", PS_LOG_WARN, "measure galaxy shapes : %f sec for %ld objects\n", psTimerMark ("psphot.galaxy"), sources->n); return true; } @@ -152,12 +193,12 @@ bool psphotGalaxyShape_Threaded (psThreadJob *job) { - pmReadout *readout = job->args->data[0]; - psArray *sources = job->args->data[2]; - 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); + pmReadout *readout = job->args->data[0]; + 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); for (int i = 0; i < sources->n; i++) { - - // XXX what do i need for the kron flux (given a supplied kron radius?) pmSource *source = sources->data[i]; @@ -184,10 +225,10 @@ // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS - bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); psAssert (source->pixels, "WTF?"); // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr // do the following for a set of shapes (Ex,Ey) - psphotGalaxyShapeSource (source, maskVal); + psphotGalaxyShapeGrid (source, fitOptions, maskVal, psfSize); psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); @@ -198,45 +239,51 @@ } -bool psphotGalaxyShapeGrid (pmSource *source, psImageMaskType maskVal, float psfSize) { - - pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize); - - // I have some source guess (e0, e1, e2) - psEllipsePol guessPol; - - // 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) - - for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) { - for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) { +bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize) { + + // allocate the model + // XXX need to figure out how to get the model into the program + pmModelType modelType = pmModelClassGetType ("PS_MODEL_EXP"); + + pmModel *model = pmModelAlloc(modelType); + if (!model) { + return NULL; + } + + pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + + psF32 *PAR = pcm->modelConv->params->data.F32; + + // I have some source guess (e0, e1, e2) + psEllipsePol guessPol; + + // 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) + + for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) { + for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) { - psEllipsePol testPol = guessPol; - testPol.e1 += de1; - testPol.e2 += de2; - - psEllipseShape shape; - psEllipsePolToShape (&shape, testPol); - - PAR[PM_PAR_SXX] = shape.sxx; - PAR[PM_PAR_SXY] = shape.sxy; - PAR[PM_PAR_SYY] = shape.syy; - - // where do the results go?? - psphotGalaxyShapeSource (pcm, source, maskVal); - } - } - - return true; + psEllipsePol testPol = guessPol; + testPol.e1 += de1; + testPol.e2 += de2; + + psEllipseAxes axes = psEllipsePolToAxes (testPol, 0.1); + + pmPSF_AxesToModel (PAR, axes, modelType); + + // where do the results go?? + psphotGalaxyShapeSource (pcm, source, maskVal, psfSize); + } + } + + return true; } // fit the given model to the source and find chisq & normalization // XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta? -bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal) { +bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize) { PS_ASSERT_PTR_NON_NULL(source, false); PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_ASSERT_PTR_NON_NULL(source->pixels, false); - - psF32 *PAR = pcm->modelConv->params->data.F32; // generate the modelFlux @@ -249,26 +296,26 @@ 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; + 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; } - // 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; - } } @@ -276,4 +323,6 @@ float Chisq = YY - 2 * Io * YM + Io * Io * MM; - return true; -} + fprintf (stderr, "Io: %f, Chisq: %f\n", Io, Chisq); + + return true; +} Index: branches/eam_branches/ipp-20130904/psphot/src/psphotKronFlux.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotKronFlux.c (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotKronFlux.c (revision 36148) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) +bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -35,8 +35,5 @@ psAssert (sources, "missing sources?"); - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - // psAssert (psf, "missing psf?"); - - if (!psphotKronFluxReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) { + if (!psphotKronFluxReadout (config, recipe, view, filerule, readout, sources)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; @@ -46,5 +43,5 @@ } -bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) { +bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources) { bool status = false; @@ -61,14 +58,4 @@ if (!status) { nThreads = 0; - } - - bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT"); - if (!status) { - KRON_APPLY_WEIGHT = true; - } - - bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW"); - if (!status) { - KRON_APPLY_WINDOW = false; } @@ -105,5 +92,5 @@ // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE"); + psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_FLUX"); psArrayAdd(job->args, 1, readout); @@ -111,5 +98,4 @@ PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, pass, PS_TYPE_S32); // set this to 0 to run without threading @@ -153,11 +139,9 @@ pmReadout *readout = job->args->data[0]; - psArray *sources = job->args->data[2]; - 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); + psArray *sources = job->args->data[1]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); for (int i = 0; i < sources->n; i++) { - - // XXX what do i need for the kron flux (given a supplied kron radius?) pmSource *source = sources->data[i]; @@ -180,9 +164,9 @@ } - // is the right?? - float windowRadius = 6.0*source->moments->Mrf ; + // is this right?? + float windowRadius = 2.5*source->moments->Mrf ; // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS - bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); psAssert (source->pixels, "WTF?"); @@ -205,4 +189,5 @@ // the peak position is less accurate but less subject to extreme deviations + // XXX do this or just use the peak position? float dX = source->moments->Mx - source->peak->xf; float dY = source->moments->My - source->peak->yf; @@ -211,7 +196,12 @@ float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + // Calculate the Kron magnitude (make this block optional?) // XXX set the aperture here - float radKron = 2.5*Mrf; + float radKron = 2.5*source->moments->Mrf; float radKron2 = radKron*radKron; @@ -222,7 +212,7 @@ // set vPix to the source pixels (it may have been set to the // smoothed image above) - vPix = source->pixels->data.F32; - vWgt = source->variance->data.F32; - vMsk = source->mask->data.F32; + psF32 **vPix = source->pixels->data.F32; + psF32 **vWgt = source->variance->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; for (psS32 row = 0; row < source->pixels->numRows ; row++) { Index: branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c (revision 36148) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) +bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; @@ -35,8 +35,5 @@ psAssert (sources, "missing sources?"); - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - // psAssert (psf, "missing psf?"); - - if (!psphotPetroFluxReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) { + if (!psphotPetroFluxReadout (config, recipe, view, filerule, readout, sources)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; @@ -46,14 +43,14 @@ } -bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) { +bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources) { bool status = false; if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping kron fluxes"); + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping petrosian fluxes"); return true; } - psTimerStart ("psphot.kron"); + psTimerStart ("psphot.petro"); // determine the number of allowed threads @@ -61,14 +58,4 @@ if (!status) { nThreads = 0; - } - - bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT"); - if (!status) { - KRON_APPLY_WEIGHT = true; - } - - bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW"); - if (!status) { - KRON_APPLY_WINDOW = false; } @@ -105,5 +92,5 @@ // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE"); + psThreadJob *job = psThreadJobAlloc ("PSPHOT_PETRO_FLUX"); psArrayAdd(job->args, 1, readout); @@ -111,5 +98,4 @@ PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, pass, PS_TYPE_S32); // set this to 0 to run without threading @@ -146,5 +132,5 @@ psFree (cellGroups); - psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); + psLogMsg ("psphot.petro", PS_LOG_WARN, "measure petro fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.petro"), sources->n); return true; } @@ -153,11 +139,9 @@ pmReadout *readout = job->args->data[0]; - psArray *sources = job->args->data[2]; - 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); + psArray *sources = job->args->data[1]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); for (int i = 0; i < sources->n; i++) { - - // XXX what do i need for the kron flux (given a supplied kron radius?) pmSource *source = sources->data[i]; @@ -169,5 +153,5 @@ if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; - // XXX where are we storing the supplied kron radius? + // XXX where are we storing the supplied petro radius? if (!isfinite(source->moments->Mrf)) continue; @@ -180,9 +164,9 @@ } - // is the right?? - float windowRadius = 6.0*source->moments->Mrf ; + // is the right?? XXX need the correct location for Rpet + float windowRadius = 2.0*source->moments->Mrf ; // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS - bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); psAssert (source->pixels, "WTF?"); @@ -211,7 +195,12 @@ float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + // Calculate the Petro magnitude (make this block optional?) // XXX set the aperture here - float radPetro = 2.5*Mrf; + float radPetro = 2.0*source->moments->Mrf; float radPetro2 = radPetro*radPetro; @@ -222,7 +211,7 @@ // set vPix to the source pixels (it may have been set to the // smoothed image above) - vPix = source->pixels->data.F32; - vWgt = source->variance->data.F32; - vMsk = source->mask->data.F32; + psF32 **vPix = source->pixels->data.F32; + psF32 **vWgt = source->variance->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -253,7 +242,9 @@ // return the flux and error to parameters? - source->moments->PetroFlux = Sum; - source->moments->PetroFluxErr = sqrt(Var); - - return true; -} + // source->moments->PetroFlux = Sum; + // source->moments->PetroFluxErr = sqrt(Var); + + fprintf (stderr, "petro flux: %f +/- %f\n", Sum, sqrt(Var)); + + return true; +} Index: branches/eam_branches/ipp-20130904/psphot/src/psphotSetThreads.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotSetThreads.c (revision 36146) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotSetThreads.c (revision 36148) @@ -47,4 +47,19 @@ psFree(task); + task = psThreadTaskAlloc("PSPHOT_KRON_FLUX", 4); + task->function = &psphotKronFlux_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_PETRO_FLUX", 4); + task->function = &psphotPetroFlux_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPE", 6); + task->function = &psphotGalaxyShape_Threaded; + psThreadTaskAdd(task); + psFree(task); + task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10); task->function = &psphotBlendFit_Threaded;