Index: branches/eam_branches/ipp-20140904/psphot/src/Makefile.am =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/Makefile.am (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/Makefile.am (revision 37804) @@ -260,5 +260,6 @@ psphotSetNFrames.c \ psphotSourceMemory.c \ - psphotChipParams.c + psphotChipParams.c \ + psphotGalaxyParams.c # not currently used Index: branches/eam_branches/ipp-20140904/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphot.h (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphot.h (revision 37804) @@ -138,9 +138,13 @@ bool psphotExtendedSourceFits_Threaded (psThreadJob *job); +bool psphotGalaxyParams (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotGalaxyParamsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotGalaxyParams_Threaded (psThreadJob *job); + bool psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotMagnitudes (pmConfig *config, const pmFPAview *view, const char *filerule); -bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf); +bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf, int index); bool psphotMagnitudes_Threaded (psThreadJob *job); @@ -168,4 +172,5 @@ // in psphotChoosePSF.c: bool psphotPSFstats (pmReadout *readout, pmPSF *psf); +bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf); bool psphotMomentsStats (pmReadout *readout, psArray *sources); @@ -232,6 +237,6 @@ int psphotSaveImage (psMetadata *header, psImage *image, char *filename); bool psphotDumpConfig (pmConfig *config); -pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view); -pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view); +pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view, int index); +pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view, int index); // PSF / DBL / EXT evaluation functions Index: branches/eam_branches/ipp-20140904/psphot/src/psphotChooseAnalysisOptions.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotChooseAnalysisOptions.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotChooseAnalysisOptions.c (revision 37804) @@ -606,4 +606,10 @@ psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld objects: %f sec\n", objects->n, psTimerMark ("psphot.options")); + psFree(exptime); + psFree(zeropt); + psFree(chips); + psFree(inputFilters); + psFree(petroFluxLim); + psFree(extFitFluxLim); psFree(toGal); Index: branches/eam_branches/ipp-20140904/psphot/src/psphotChoosePSF.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotChoosePSF.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotChoosePSF.c (revision 37804) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf); // generate a PSF model for inputs without PSF models already loaded @@ -388,5 +389,7 @@ psFree (models); - if (!psphotPSFstats (readout, psf)) { + // if (!psphotPSFstats (readout, psf)) { + + if (!psphotPSFstatsSources (readout, sources, psf)) { psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); psFree(options); @@ -466,13 +469,13 @@ if (modelPSF->params->n > 7) { - psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]); + psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]); } if (modelPSF->params->n > 8) { - psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]); + psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]); } psFree (modelPSF); - // float fwhmtest = pmPSFtoFWHM(psf, xc, yc); - // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); + float fwhmtest = pmPSFtoFWHM(psf, xc, yc); + fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); } } @@ -522,19 +525,172 @@ if (psfExtra1->n) { - if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1"); + if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev); + } + + if (psfExtra2->n) { + if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev); + } + + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars); + + // psLogMsg ("psphot", PS_LOG_DETAIL, "PSF angle: %f, nstars: %d", axes.theta, psf->nPSFstars); + + char *psfModelName = pmModelClassGetName(psf->type); + psMetadataAddStr(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", psfModelName); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK", PS_META_REPLACE, "Valid PSF Model?", true); + + int nParams = pmModelClassParameterCount(psf->type); + psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "PSF_NPAR", PS_META_REPLACE, "Number of PSF parameters", nParams); + + psFree (fwhmMajor); + psFree (fwhmMinor); + psFree (psfExtra1); + psFree (psfExtra2); + psFree (stats); + return true; + +escape: + psFree (fwhmMajor); + psFree (fwhmMinor); + psFree (psfExtra1); + psFree (psfExtra2); + psFree (stats); + return false; +} + +// measure average parameters of the PSF model +bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf) { + + psEllipseShape shape; + psEllipseAxes axes; + + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(psf, false); + + psImage *image = readout->image; + PS_ASSERT_PTR_NON_NULL(image, false); + + // XXX a minor hack: measure the values for a grid of points across the image, determine mean, UQ, LQ, etc: + psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *psfExtra1 = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *psfExtra2 = psVectorAllocEmpty (100, PS_DATA_F32); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + float xc = source->peak->xf; + float yc = source->peak->yf; + + // create modelPSF from this model + pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0); + if (!modelPSF) { + fprintf (stderr, "?"); + continue; + } + + // get the model full-width at half-max + float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5); + + // XXX make sure this is consistent with the re-definition of PM_PAR_SXX + shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; + shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; + shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; + axes = psEllipseShapeToAxes (shape, 20.0); + + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) { + fprintf (stderr, "!"); + psFree (modelPSF); + continue; + } + psVectorAppend (fwhmMajor, FWHM_MAJOR); + psVectorAppend (fwhmMinor, FWHM_MINOR); + + if (modelPSF->params->n > 7) { + psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]); + } + if (modelPSF->params->n > 8) { + psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]); + } + psFree (modelPSF); + + float fwhmtest = pmPSFtoFWHM(psf, xc, yc); + fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE); + if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR"); goto escape; - } - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean); - psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev); + } + + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ", PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG", PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ", PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ", PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ); + + float fwhmMaj = stats->sampleMean; // FWHM on major axis + + if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN", PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG", PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ", PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ", PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ); + + float fwhmMin = stats->sampleMean; // FWHM on minor axis + + if (readout->parent) { + + // we now have 2 possible measurements of the seeing : the PSF model based version and + // the Moments based value we need to define a definitive "CHIP.SEEING" value. This + // value is used by the PSF-matching programs (ppSub and ppStack). The moments-based + // value is probably more representative of the value needed by those tools (it is also + // more stable?) + bool status = false; + float fwhmMajorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1"); + float fwhmMinorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1"); + + pmChip *chip = readout->parent->parent; // Parent chip + psAssert(chip, "Cell should be attached to a chip."); + psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip + item->data.F32 = 0.5 * (fwhmMajorMoments + fwhmMinorMoments); + + psLogMsg ("psphot", PS_LOG_DETAIL, "fwhm (psf): %f,%f (moments): %f,%f", fwhmMaj, fwhmMin, fwhmMajorMoments, fwhmMinorMoments); + } + + if (psfExtra1->n) { + if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev); } if (psfExtra2->n) { - if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2"); - goto escape; - } - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean); - psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev); + if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev); } Index: branches/eam_branches/ipp-20140904/psphot/src/psphotCleanup.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotCleanup.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotCleanup.c (revision 37804) @@ -25,5 +25,5 @@ // don't display memory leaks unless trace is on. // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot"); - fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot"); + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks2 (0, NULL, stdout, false, 500), "psphot"); #endif return; Index: branches/eam_branches/ipp-20140904/psphot/src/psphotEfficiency.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotEfficiency.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotEfficiency.c (revision 37804) @@ -502,5 +502,5 @@ // measure the magnitudes and fluxes for the sources - if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) { + if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf, index)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); psFree(fakeSources); Index: branches/eam_branches/ipp-20140904/psphot/src/psphotExtendedSourceFits.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotExtendedSourceFits.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotExtendedSourceFits.c (revision 37804) @@ -629,5 +629,5 @@ pmSourceCacheModel (source, maskVal); -# if (PS_TRACE_ON) +# if (0 && PS_TRACE_ON) pmModel *model = source->modelFits->data[0]; int flags = 0xffffffff; Index: branches/eam_branches/ipp-20140904/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotFitSourcesLinear.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotFitSourcesLinear.c (revision 37804) @@ -263,4 +263,10 @@ } } + if (modelSum == 0 || maskedSum == 0) { + // no pixels at all drop the modelPSF (probably only an issue for external sources + // which are created from a fake peak (full force)) + psFree(source->modelPSF); + continue; + } if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit) if (modelSum < 0.8) { Index: branches/eam_branches/ipp-20140904/psphot/src/psphotFullForceReadout.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotFullForceReadout.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotFullForceReadout.c (revision 37804) @@ -117,4 +117,10 @@ psphotGalaxyShape (config, view, filerule); + // measure aperture photometry corrections + if (!psphotApResid (config, view, filerule)) { + psLogMsg ("psphot", 3, "failed on psphotApResid"); + return psphotReadoutCleanup (config, view, filerule); + } + // calculate source magnitudes (psf mag and ap mag) psphotMagnitudes(config, view, filerule); Index: branches/eam_branches/ipp-20140904/psphot/src/psphotGalaxyParams.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotGalaxyParams.c (revision 37804) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotGalaxyParams.c (revision 37804) @@ -0,0 +1,637 @@ +# include "psphotInternal.h" + +#define USE_FOOTPRINTS 1 + +// measure some properties on the exended objects + +bool psphotGalaxyParams (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Galaxy Parameters ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // add this + bool doGalaxyParams = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_GALAXY_PARAMS"); + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doExtModels = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS"); + + // All of these must be true to measure galaxy parameters + if (!doGalaxyParams || !doPetrosian || !doAnnuli || !doExtModels) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping Galaxy parameter source measurements\n"); + return true; + } + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotGalaxyParamsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on galaxy parameter measurements for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +#ifdef notdef +/*** for the moment, this test code : it is not thread safe ***/ +int Nall = 0; +int Nskip1 = 0; +int Nskip2 = 0; +int Nskip3 = 0; +int Nskip4 = 0; +int Nskip5 = 0; +int Nskip6 = 0; + +# define SKIP(VALUE) { VALUE++; continue; } + +#endif //notdef + +#define SKIP(VALUE) {continue;} + +// aperture-like measurements for extended sources +bool psphotGalaxyParamsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + int Next = 0; + int Npetro = 0; + int Nannuli = 0; + + psTimerStart ("psphot.galaxy"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // get the sky noise from the background analysis; if this is missing, get the user-supplied value + float skynoise = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV"); + if (!status) { + skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + psWarning ("failed to get sky noise level from background analysis; defaulting to user supplied value of %f\n", skynoise); + } + + float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio + if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) { + fitNsigmaConv = 5.0; + } + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + + psImage *footprintImage = NULL; +#ifdef USE_FOOTPRINTS + footprintImage = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32); + psImageInit(footprintImage, 0); + pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, true); +#endif + + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0); + *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined"); + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_PARAMS"); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, AnalysisRegion); + psArrayAdd(job->args, 1, recipe); + + PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, fitNsigmaConv, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_F32); + psArrayAdd(job->args, 1, footprintImage); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli + + +// set this to 0 to run without threading +# if (0) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } +# else + if (!psphotGalaxyParams_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + psScalar *scalar = NULL; + scalar = job->args->data[8]; + Next += scalar->data.S32; + scalar = job->args->data[9]; + Npetro += scalar->data.S32; + scalar = job->args->data[10]; + Nannuli += scalar->data.S32; + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[8]; + Next += scalar->data.S32; + scalar = job->args->data[9]; + Npetro += scalar->data.S32; + scalar = job->args->data[10]; + Nannuli += scalar->data.S32; + } + psFree(job); + } + } + psFree (cellGroups); + psFree(AnalysisRegion); + psFree(footprintImage); + + psLogMsg ("psphot", PS_LOG_WARN, "galaxy parameter analysis: %f sec for %d objects\n", psTimerMark ("psphot.galaxy"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); + psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); + +// # if (PS_TRACE_ON) +# if (0) + fprintf (stderr, "ext analysis skipped @ 1 : %d\n", Nskip1); + fprintf (stderr, "ext analysis skipped @ 2 : %d\n", Nskip2); + fprintf (stderr, "ext analysis skipped @ 3 : %d\n", Nskip3); + fprintf (stderr, "ext analysis skipped @ 4 : %d\n", Nskip4); + fprintf (stderr, "ext analysis skipped @ 5 : %d\n", Nskip5); + fprintf (stderr, "ext analysis skipped @ 6 : %d\n", Nskip6); +#endif + + // psphotRadialProfileShowSkips (); + + // psphotVisualShowResidualImage (readout, false); + + // psphotVisualShowPetrosians (sources); + + return true; +} + +bool psphotGalaxyParams_Threaded (psThreadJob *job) { + + // bool status; + + int Next = 0; + int Npetro = 0; + int Nannuli = 0; + int Ngood = 0; + + // arguments: readout, sources, models, region, psfSize, maskVal, markVal + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + psRegion *region = job->args->data[2]; + psMetadata *recipe = job->args->data[3]; + +#ifndef USE_FOOTPRINTS + float skynoise = PS_SCALAR_VALUE(job->args->data[4],F32); +#endif + float fitNsigmaConv = PS_SCALAR_VALUE(job->args->data[5],F32); + float psfSize = PS_SCALAR_VALUE(job->args->data[6],F32); + psImage *footprintImage = job->args->data[7]; + + bool status; +// bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // don't try and use bad models + pmModelStatus badModel = PM_MODEL_STATUS_NONE; + badModel |= PM_MODEL_STATUS_NONCONVERGE; + badModel |= PM_MODEL_STATUS_BADARGS; + badModel |= PM_MODEL_STATUS_OFFIMAGE; + badModel |= PM_MODEL_STATUS_NAN_CHISQ; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID; + badModel |= PM_MODEL_PCM_FAIL_GUESS; + + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + + psVector *bgPixels = NULL; + psVector *bgPixelsFlip = NULL; +#define dontMAKE_MARK_IMAGE 1 +#ifdef MAKE_MARK_IMAGE + psImage *markImage = NULL; +#endif + + #define MARK_SOURCE 0 + #define MARK_MASK 1 + #define MARK_BG 2 + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + +#define dontDUMP_IMAGES 1 +#ifdef DUMP_IMAGES + char *filebase = "/data/ippc18.0/home/bills/debug/psphot/images"; + char filename[80]; + psMetadata *dummyHeader = psMetadataAlloc(); +#endif + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // Nall ++; + + // rules for measuring petrosian parameters for specific objects are set in + // psphotChooseAnalysisOptions.c + if (!(source->tmpFlags & PM_SOURCE_TMPF_PETRO)) SKIP (Nskip1); + + // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?) + if (source->peak->x < region->x0) SKIP (Nskip2); + if (source->peak->y < region->y0) SKIP (Nskip3); + if (source->peak->x > region->x1) SKIP (Nskip4); + if (source->peak->y > region->y1) SKIP (Nskip5); + + // Should we skip objects that do not have an extended model? + if (!source->extpars) SKIP (Nskip6); + if (!source->modelFits) SKIP (Nskip7); + + // Flux cut to reduce sources to those that match SDSS + // if (source->moments->KronFlux < 2e5) continue; + + pmModel *sersicModel = NULL; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + if (model->type == sersicModelType) { + sersicModel = model; + break; + } + } + if (!sersicModel) SKIP (Nskip8); + if (sersicModel->flags & badModel) SKIP (Nskip8); + Next ++; + + // if sersicModel is not modelEXT make it so + pmSourceMode saveMode = source->mode; + pmSourceType saveType = source->type; + bool isExtended = saveType == PM_SOURCE_TYPE_EXTENDED; + pmModel *extModel = source->modelEXT; + + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + if (!isExtended || extModel != sersicModel) { + // XXX: I'm not doing this right. For now just look at objects for which + source->modelEXT = psMemIncrRefCounter (sersicModel); + source->type = PM_SOURCE_TYPE_EXTENDED; + source->mode |= PM_SOURCE_MODE_EXTMODEL; + // do we need this? + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; +#ifdef notdef +#endif + } + + + psF32 *PAR = sersicModel->params->data.F32; + psF32 Xo = PAR[PM_PAR_XPOS]; + psF32 Yo = PAR[PM_PAR_YPOS]; + psF32 Io = PAR[PM_PAR_I0]; + // For now: set the half light radius to the major axis of the model + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff); + source->extpars->ghalfLightRadius = axes.major; + // and we will integrate out to 2 * half light radius + psF32 R = 2 * axes.major; + psF32 R2 = PS_SQR(R); + + float radius = source->peak->xf - source->pixels->col0; + radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + // Make image at least 10 * half light radius + radius = PS_MAX (radius, 5*axes.major); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 2*radius); + + // RedefinePixels frees modelFlux + if (source->modelEXT->isPCM) { + pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv); + } else { + pmSourceCacheModel (source, maskVal); + } + +#ifdef DUMP_IMAGES + sprintf(filename, "%s/%04d-pix.fits", filebase, source->id); + + psphotSaveImage(dummyHeader, source->pixels, filename); +#endif + + + // Create the markImage and accumlate the Flux + + // 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 + + psF32 **vPix = source->pixels->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + bgPixels = psVectorRecycle(bgPixels, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32); + bgPixels->n = 0; + bgPixelsFlip = psVectorRecycle(bgPixelsFlip, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32); + bgPixelsFlip->n = 0; + +#ifdef MAKE_MARK_IMAGE + // The markImage may no longer be needed + markImage = psImageRecycle(markImage, source->pixels->numCols, source->pixels->numRows, PS_TYPE_U8); + + // assume all pixels are masked to start + psImageInit(markImage, MARK_MASK); + psU8 **vMark = markImage->data.U8; + #define SET_MARK(_val) vMark[row][col] = _val +#else + #define SET_MARK(_val) +#endif + +#ifdef USE_FOOTPRINTS + psS32 **vFootprint = footprintImage->data.S32; +#else + // XXX use recipe value instead of hard coded 1.5 + psF32 threshold = skynoise * 3; +#endif + + psF32 sumI = 0; + psS32 row0 = source->pixels->row0; + psS32 col0 = source->pixels->col0; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + psF32 yDiff = row - yCM; + // y coordinate of mirror pixel + int yFlip = yCM - yDiff; + if (yFlip < 0) continue; + if (yFlip >= source->pixels->numRows) continue; + psS32 yImage = row + row0; + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + psF32 pixel = vPix[row][col]; + if (isnan(pixel)) continue; + // valid pixel + + psS32 xImage = col + col0; + + psF32 xDiff = col - xCM; + + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 <= R2) { + sumI += pixel; + } + +#ifdef USE_FOOTPRINTS + if (vFootprint[yImage][xImage]) { +#else + if (fabs(pixel) >= threshold) { +#endif + SET_MARK(MARK_SOURCE); + continue; + } + // we have a background pixel. + // Check whether it's mirror is a background pixel as well. + + // x coordinate of mirror pixel + int xFlip = xCM - xDiff; + if (xFlip < 0) continue; + if (xFlip >= source->pixels->numCols) continue; + // check mask and value for mirror pixel + if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue; + psF32 mirrorPixel = vPix[yFlip][xFlip]; + if (isnan(mirrorPixel)) continue; + +#ifdef USE_FOOTPRINTS + if (vFootprint[yFlip + row0][xFlip + col0] == 0) { +#else + if (fabs(mirrorPixel) < threshold) { +#endif + // This is a good background pixel with a good mirror Save the values + SET_MARK(MARK_BG); + psVectorAppend(bgPixels, pixel); + psVectorAppend(bgPixelsFlip, mirrorPixel); + } + } + } +#ifdef DUMP_IMAGES + sprintf(filename, "%s/%04d-mrk.fits", filebase, source->id); + + psphotSaveImage(dummyHeader, markImage, filename); +#endif + + // subtract model to create the residual image + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + // Hmm something has gone wrong + fprintf(stderr, "failed to subtract sersic model for source %d\n", source->id); + goto restoreModel; + } +#ifdef DUMP_IMAGES + sprintf(filename, "%s/%04d-res.fits", filebase, source->id); + + psphotSaveImage(dummyHeader, source->pixels, filename); +#endif + + // don't know if these pointers could change in the background subtraction so grab them again + vPix = source->pixels->data.F32; + vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + psF32 **vVar = source->variance->data.F32; + psF32 **vModel = source->modelFlux->data.F32; + // Accumulate Rt and Ra + psF32 sumRt = 0; + psF32 sumRa = 0; + int numPixels = 0; + psF32 sumB = 0; + psF32 sumModel = 0; + int numB = 0; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + psF32 yDiff = row - yCM; + if (fabs(yDiff) > R) continue; + // y coordinate of mirror pixel + int yFlip = yCM - yDiff; + if (yFlip < 0) continue; + if (yFlip >= source->pixels->numRows) continue; + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + psF32 pixel = vPix[row][col]; + if (isnan(pixel)) continue; + + psF32 xDiff = col - xCM; + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > R2) continue; + + // Measurement of bumpiness parameter excludes the pixels + // within 2 pixels of the center + if (r2 > 4) { + psF32 variance = vVar[row][col]; + if (isfinite(variance)) { + psF32 modelPixel = vModel[row][col]; + if (!isfinite(modelPixel)) { + // XXX: Can this happen in a good model? + fprintf(stderr, "non-finite modelPixel %4d (%4d, %4d)\n", source->id, col, row); + } + numB++; + sumB += PS_SQR(pixel) - variance; + sumModel += Io * modelPixel; + } + } + + // x coordinate of mirror pixel + int xFlip = xCM - xDiff; + if (xFlip < 0) continue; + if (xFlip >= source->pixels->numCols) continue; + // check mask and value for mirror pixel + if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue; + psF32 mirrorPixel = vPix[yFlip][xFlip]; + if (isnan(mirrorPixel)) continue; + + sumRt += fabs(pixel + mirrorPixel); + sumRa += fabs(pixel - mirrorPixel); + ++numPixels; + } + } + + if (numB > 0) { + source->extpars->gbumpy = sqrt(sumB/numB) / (sumModel / numB); + } + + if (bgPixels->n < numPixels) { + fprintf(stderr, "could not find enough background pixels for source %d found %ld need %d\n", source->id, bgPixels->n, numPixels); + goto restoreModel; + } + Ngood++; + + // and calculate Bt and Ba + psF32 sumBt = 0; + psF32 sumBa = 0; + psF32 *vB = bgPixels->data.F32; + psF32 *vBFlip = bgPixelsFlip->data.F32; + long length = bgPixels->n; + // Randomly select numPixels background pixels + for (int n = 0; n < numPixels; n++) { + double rand = psRandomUniform(rng); + long index = floor(rand * length); + if (index >= length) { + fprintf(stderr, "random index %ld is larger than vector length %ld\n", + index, length); + continue; + } + sumBt += fabs(vB[index] + vBFlip[index]); + sumBa += fabs(vB[index] - vBFlip[index]); + } + + // finally save the results + source->extpars->gRT = 0.5 * (sumRt - sumBt)/sumI; + source->extpars->gRA = 0.5 * (sumRa - sumBa)/sumI; + source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA; + +#ifdef PRINT_STUFF + if (Ngood % 40 == 1) { + printf(" id sumI rHL rad sumI gRT gRA sumRt sumBt numPixels numBGPixels\n"); + } + printf("%5d %5.2e %3.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n", + source->id, sumI, source->extpars->ghalfLightRadius, radius, source->extpars->gRT, source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n); +#endif + + // Now put things back the way that they were. +restoreModel: + + if ((extModel != sersicModel) || !isExtended) { + // add the sersic model back in + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // replace the extended model + if (extModel != sersicModel) { + psMemDecrRefCounter (sersicModel); + source->modelEXT = extModel; + } + source->mode = saveMode; + source->type = saveType; + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + } + +#ifdef MAKE_MARK_IMAGE + psFree(markImage); +#endif +#ifdef DUMP_IMAGES + psFree(dummyHeader); + // psphotSaveImage raises an error + psErrorClear(); +#endif + psFree(bgPixels) + psFree(bgPixelsFlip) + psFree(rng); + + psScalar *scalar = NULL; + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + scalar = job->args->data[8]; + scalar->data.S32 = Next; + + scalar = job->args->data[9]; + scalar->data.S32 = Npetro; + + scalar = job->args->data[10]; + scalar->data.S32 = Nannuli; + + return true; +} + Index: branches/eam_branches/ipp-20140904/psphot/src/psphotMagnitudes.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotMagnitudes.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotMagnitudes.c (revision 37804) @@ -38,5 +38,5 @@ psAssert (psf, "missing psf?"); - if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) { + if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; @@ -46,5 +46,5 @@ } -bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) { +bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf, int index) { bool status = false; @@ -65,6 +65,6 @@ // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN - pmReadout *backModel = psphotSelectBackground (config, view); - pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); + pmReadout *backModel = psphotSelectBackground (config, view, index); + pmReadout *backStdev = psphotSelectBackgroundStdev (config, view, index); // bit-masks to test for good/bad pixels Index: branches/eam_branches/ipp-20140904/psphot/src/psphotOutput.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotOutput.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotOutput.c (revision 37804) @@ -29,10 +29,7 @@ } -pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view) { - - bool status; - - - pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL")); +pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view, int index) { + + pmFPAfile *file = pmFPAfileSelectSingle (config->files, psphotGetFilerule("PSPHOT.BACKMDL"), index); if (!file) return NULL; @@ -41,9 +38,7 @@ } -pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view) { - - bool status; - - pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV")); +pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view, int index) { + + pmFPAfile *file = pmFPAfileSelectSingle (config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"), index); if (!file) return NULL; Index: branches/eam_branches/ipp-20140904/psphot/src/psphotPetrosianAnalysis.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotPetrosianAnalysis.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotPetrosianAnalysis.c (revision 37804) @@ -12,11 +12,4 @@ // XXX temporary user-supplied systematic sky noise measurement (derive from background model) float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); - -# if (0) - // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN - // XXX use this to set skynoise - pmReadout *backModel = psphotSelectBackground (config, view); - pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); -# endif // S/N limit to perform full non-linear fits Index: branches/eam_branches/ipp-20140904/psphot/src/psphotReadout.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotReadout.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotReadout.c (revision 37804) @@ -325,4 +325,6 @@ psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources) psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources) + // measure some parameters for galaxy science + psphotGalaxyParams (config, view, filerule); psphotRadialApertures(config, view, filerule, 0); Index: branches/eam_branches/ipp-20140904/psphot/src/psphotReadoutCleanup.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotReadoutCleanup.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotReadoutCleanup.c (revision 37804) @@ -57,5 +57,5 @@ // use the psf-model to measure FWHM stats if (psf) { - if (!psphotPSFstats (readout, psf)) { + if (!psphotPSFstatsSources (readout, sources, psf)) { psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters"); return false; Index: branches/eam_branches/ipp-20140904/psphot/src/psphotSetNFrames.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotSetNFrames.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotSetNFrames.c (revision 37804) @@ -45,8 +45,14 @@ pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - psAssert (detections->allSources, "missing sources?"); + if (!detections) { + psLogMsg ("psphot", PS_LOG_INFO, "No detections image input %d", index); + return true; + } psArray *sources = detections->allSources; + if (!sources) { + psLogMsg ("psphot", PS_LOG_INFO, "No sources for input %d", index); + return true; + } int col0 = expnum->col0; Index: branches/eam_branches/ipp-20140904/psphot/src/psphotSetThreads.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotSetThreads.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotSetThreads.c (revision 37804) @@ -92,4 +92,9 @@ psFree(task); + task = psThreadTaskAlloc("PSPHOT_GALAXY_PARAMS", 11); + task->function = &psphotGalaxyParams_Threaded; + psThreadTaskAdd(task); + psFree(task); + return true; } Index: branches/eam_branches/ipp-20140904/psphot/src/psphotSourceStats.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotSourceStats.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotSourceStats.c (revision 37804) @@ -185,18 +185,18 @@ psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS"); - psArrayAdd(job->args, 1, cells->data[j]); // sources - - PS_ARRAY_ADD_SCALAR(job->args, INNER, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, MIN_SN, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, SIGMA, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, 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); - - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint + psArrayAdd(job->args, 1, cells->data[j]); // 0 : sources + + PS_ARRAY_ADD_SCALAR(job->args, INNER, PS_TYPE_F32); // 1 + PS_ARRAY_ADD_SCALAR(job->args, MIN_SN, PS_TYPE_F32); // 2 + PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); // 3 + PS_ARRAY_ADD_SCALAR(job->args, SIGMA, PS_TYPE_F32); // 4 + PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); // 5 + + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); // 6 + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); // 7 + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // 8 : this is used as a return value for Nmoments + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // 9 : this is used as a return value for Nfail + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // 10 : this is used as a return value for Nfaint if (!psThreadJobAddPending(job)) { @@ -316,4 +316,5 @@ // XXX: this must match the above + psAbort ("this code is broken: must match psphotSourceStatsReadout above"); psArrayAdd(job->args, 1, cells->data[j]); // sources psArrayAdd(job->args, 1, recipe); @@ -472,11 +473,11 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[7]; + scalar = job->args->data[8]; scalar->data.S32 = Nmoments; - scalar = job->args->data[8]; + scalar = job->args->data[9]; scalar->data.S32 = Nfail; - scalar = job->args->data[9]; + scalar = job->args->data[10]; scalar->data.S32 = Nfaint; Index: branches/eam_branches/ipp-20140904/psphot/src/psphotStackParseCamera.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotStackParseCamera.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotStackParseCamera.c (revision 37804) @@ -262,4 +262,5 @@ } +#ifdef MAKE_CHISQ_IMAGE // generate an pmFPAimage for the chisqImage { @@ -293,4 +294,7 @@ chisqVariance->save = saveChisq; } +#else + (void) saveChisq; +#endif // MAKE_CHISQ_IMAGE psTrace("psphot", 1, "Done with psphotStackParseCamera...\n"); Index: branches/eam_branches/ipp-20140904/psphot/src/psphotStackReadout.c =================================================================== --- branches/eam_branches/ipp-20140904/psphot/src/psphotStackReadout.c (revision 37363) +++ branches/eam_branches/ipp-20140904/psphot/src/psphotStackReadout.c (revision 37804) @@ -341,4 +341,7 @@ // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N) psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources) + + // measure some parameters for galxy science + psphotGalaxyParams (config, view, STACK_RAW); // pass 1 (detections->allSources) // create source children for the OUT filerule (for radial aperture photometry and output) @@ -354,4 +357,7 @@ return psphotReadoutCleanup (config, view, STACK_RAW); } + // These arrays are no longer needed. The inputs' sources arrays have references to the data + psFree (objects); + psFree (objectsOut); bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES"); @@ -410,6 +416,4 @@ // measure aperture photometry corrections if (!psphotApResid (config, view, STACK_RAW)) { - psFree (objects); - psFree (objectsOut); psLogMsg ("psphot", 3, "failed on psphotApResid"); return psphotReadoutCleanup (config, view, STACK_RAW); @@ -432,5 +436,5 @@ logMemStats("final"); -#if (1) +#if (0) psphotSourceMemory(config, view, STACK_RAW); psphotSourceMemory(config, view, STACK_OUT); @@ -448,7 +452,4 @@ psphotStackRemoveChisqFromInputs(config, STACK_RAW); #endif - - psFree (objects); - psFree (objectsOut); // create the exported-metadata and free local data