Index: /branches/eam_branches/ipp-20130904/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36145) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36146) @@ -18,4 +18,5 @@ PSPHOT_SINGLE, PSPHOT_FORCED, + PSPHOT_FULL_FORCE, PSPHOT_MAKE_PSF, PSPHOT_MODEL_TEST, Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForce.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForce.c (revision 36146) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForce.c (revision 36146) @@ -0,0 +1,35 @@ +# include "psphotStandAlone.h" + +int main (int argc, char **argv) { + + psTimerStart ("complete"); + pmErrorRegister(); // register psModule's error codes/messages + psphotInit(); + + // load command-line arguments, options, and system config data + pmConfig *config = psphotFullForceArguments (argc, argv); + assert(config); + + psphotVersionPrint(); + + // load input data (config and images (signal, noise, mask) + if (!psphotParseCamera (config)) { + psErrorStackPrint(stderr, "Error setting up the camera\n"); + exit (psphotGetExitStatus()); + } + + // call psphot for each readout + if (!psphotImageLoop (config, PSPHOT_FORCED)) { + psErrorStackPrint(stderr, "Error in the psphot image loop\n"); + exit (psphotGetExitStatus()); + } + + psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete")); + + psErrorCode exit_status = psphotGetExitStatus(); + psphotCleanup (config); + exit (exit_status); +} + +// all functions which return to this level must raise one of the top-level error codes if they +// exit with an error. these error codes are used to specify the program exit status Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceArguments.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceArguments.c (revision 36146) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceArguments.c (revision 36146) @@ -0,0 +1,205 @@ +# include "psphotStandAlone.h" + +static void usage(const char *program, psMetadata *arg, pmConfig *config, int exitCode); +static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile); + +pmConfig *psphotFullForceArguments(int argc, char **argv) { + + int N; + bool status, status1, status2; + + // load config data from default locations + pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE); + if (config == NULL) { + psErrorStackPrint(stderr, "Can't read site configuration"); + exit(PS_EXIT_CONFIG_ERROR); + } + + PS_ARGUMENTS_GENERIC (psphot, config, argc, argv); + + // save the following additional recipe values based on command-line options + // these options override the PSPHOT recipe values loaded from recipe files + psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE); + + // Number of threads is handled + PS_ARGUMENTS_THREADS( psphot, config, argc, argv ); + + if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(argv[0], config, stdout); + if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(argv[0], config, stdout); + + // visual : interactive display mode + if ((N = psArgumentGet (argc, argv, "-visual"))) { + psArgumentRemove (N, &argc, argv); + pmVisualSetVisual(true); + } + + // break : used from recipe throughout psphotReadout + if ((N = psArgumentGet (argc, argv, "-break"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF + if ((N = psArgumentGet (argc, argv, "-region"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // chip selection is used to limit chips to be processed + if ((N = psArgumentGet (argc, argv, "-chip"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // if these command-line options are supplied, load the file name lists into config->arguments + // override any configuration-specified source for these files + // + pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist"); + pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist"); + + // XXX make psf model optional? + status = pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); + if (!status) { + psError(PSPHOT_ERR_ARGUMENTS, true, "No psf model is supplied (use -psf)"); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + status1 = pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC", "-src", "-srclist"); + status2 = pmConfigFileSetsMD (config->arguments, &argc, argv, "SRCTEXT", "-srctext", "-srctextlist"); + + if (!status1 && !status2) { + psError(PSPHOT_ERR_ARGUMENTS, true, "No source list is supplied (use one of -src, -srctext, -srclist, or -srctextlist)"); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + if (argc == 1) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + // the input file is a required argument; if not found, we will exit + status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); + if (!status) { + psError(PSPHOT_ERR_ARGUMENTS, false, "pmConfigFileSetsMD failed to parse arguments"); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + if (argc != 2) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Expected to see one more argument; saw %d", argc - 1); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + // output position is fixed + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + + psTrace("psphot", 1, "Done with psphotFullForceArguments...\n"); + return (config); +} + +static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile) +{ + fprintf(ofile, + "Usage: one of the following\n" + "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n" + " -variance varfile1[,varfile2,...] OutFileBaseName\n" + "\n" + "%s -list FileNameList [-masklist MaskFileNameList] \n" + " -variancelist VarFileNameList OutFileBaseName\n" + "\n" + "%s -help\n" + "\n" + "%s -version\n" + "\n" + "where:\n" + " FileNameList is a text file containing filenames, one per line\n" + " MaskFileNameList is a text file of mask filenames, one per line\n" + " VarFileNameList is a text file of variance filenames, one per line\n" + " OutFileBaseName is the 'root name' for output files\n" + "also required:\n" + " -src SrcFile1[,SrcFile2,...] or -srclist SrcFileNameList\n" + " specify sources to be measured (required)\n" + "\n" + "additional options:\n" + " -psf PsfFile1[,PsfFile2,...] or -psflist PsfFileNameList\n" + " specify PSF rather than letting %s estimate it\n" + " -chip nn[,nn,...]\n" + " select detector chips to process; default is all.\n" + " Indices correspond to zero-based offset in the FPA metadata table.\n" + " -region RegionString\n" + " specify analysis region. String is of form '[x0:x1,y0:y1]'\n" + " To use this option you must define a default in psphot.config\n" + " -visual\n" + " turns on interactive display mode\n" + " -dumpconfig CfgFileName\n" + " causes config info to be dumped to the named file.\n" + " -break NOTHING|BACKMDL|PEAKS|MOMENTS|PSFMODEL|ENSEMBLE|PASS1\n" + " choose a point at which to exit processing early\n" + " -threads n\n" + " set number of parallel threads of execution\n" + " -F OldFileRule ReplacementFileRule\n" + " change file naming rule; e.g. '-F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF'\n" + " -D name stringval\n" + " set a string-valued config parameter\n" + " -Di name intval\n" + " set an integer-valued config parameter\n" + " -Df name fval\n" + " set a float-valued config parameter\n" + " -Db name boolval\n" + " set a boolean-valued config parameter\n" + " -v, -vv, -vvv\n" + " set increasing levels of verbosity\n" + " -logfmt FormatString\n" + " set format string used for log messages\n" + " -trace Fac Lvl\n" + " set tracing for facility Fac to integer Lvl, e.g. '-trace err 10'\n" + " -trace-levels\n" + " print current trace levels\n", + program,program,program,program,program); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_SUCCESS); +} + +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config, // Configuration + int exitCode + ) +{ + fprintf(stderr, + "Usage: one of the following\n" + "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n" + " -variance varfile1[,varfile2,...] OutFileBaseName\n" + "\n" + "%s -list FileNameList [-masklist MaskFileNameList] \n" + " -variancelist VarFileNameList OutFileBaseName\n" + "also required:\n" + " -src SrcFile1[,SrcFile2,...] or -srclist SrcFileNameList\n" + " specify sources to be measured (required)\n" + "\n" + "Try '%s -help' for more options and explanation\n", + program,program,program); + if (exitCode != PS_EXIT_SUCCESS) + psErrorStackPrint(stderr, "Error reading arguments\n"); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(exitCode); +} + Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceImageLoop.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceImageLoop.c (revision 36146) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceImageLoop.c (revision 36146) @@ -0,0 +1,153 @@ +# include "psphotStandAlone.h" + +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } + +bool psphotFullForceImageLoop (pmConfig *config) { + + bool status; + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + pmFPAfile *load = psMetadataLookupPtr (&status, config->files, "PSPHOT.LOAD"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + pmHDU *lastHDU = NULL; // Last HDU updated + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + + // for psphot, we force data to be read at the chip level + while ((chip = pmFPAviewNextChip (view, load->fpa, 1)) != NULL) { + psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + + // load just the input image data (image, mask, weight) + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSPHOT.LOAD"); + pmFPAfileActivate (config->files, true, "PSPHOT.MASK"); + pmFPAfileActivate (config->files, true, "PSPHOT.VARIANCE"); + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphot."); + + // mosaic the cells of a chip into a single contiguous (trimmed) chip + if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip."); + + // Read WCS if easy. + // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell. + { + pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip + pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL); + if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) { + psWarning("Unable to read WCS astrometry from header."); + psErrorClear(); + pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL); + if (hduHigh && hduHigh != hduLow && + !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { + psWarning("Unable to read WCS astrometry from primary header."); + psErrorClear(); + } + } + } + + // try to load other supporting data (PSF, SRC, etc). + // do not re-load the following three files + pmFPAfileActivate (config->files, true, NULL); + pmFPAfileActivate (config->files, false, "PSPHOT.LOAD"); + pmFPAfileActivate (config->files, false, "PSPHOT.MASK"); + pmFPAfileActivate (config->files, false, "PSPHOT.VARIANCE"); + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphot."); + + // re-activate files so they will be closed and freed below + pmFPAfileActivate (config->files, true, NULL); + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + + // Update the header + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; + } + + // if an external mask is supplied, ensure that NAN pixels are also masked + if (readout->mask) { + psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels + if (!pmReadoutMaskInvalid(readout, maskTest, maskSat)) { + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); + psFree(view); + return false; + } + } + + // run the actual photometry analysis on this chip/cell/readout + if (!psphotFullForceReadout (config, view, "PSPHOT.INPUT")) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + } + + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } + } + + // save output which is saved at the chip level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); + } + // save output which is saved at the fpa level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); + + // fail if we failed to handle an error + if (psErrorCodeLast() != PS_ERR_NONE) psAbort ("failed to handle an error!"); + + psFree (view); + return true; +} + +// I/O files related to psphot: +// PSPHOT.INPUT : input image file(s) +// PSPHOT.RESID : residual image +// PSPHOT.OUTPUT : output object tables (object) + +// PSPHOT.BACKSUB : background subtracted image +// PSPHOT.BACKGND : background model (full-scale image?) +// PSPHOT.BACKMDL : background model (binned image?) +// PSPHOT.PSF : sample PSF images + +// PSPHOT.MASK +// PSPHOT.VARIANCE +// 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 36146) @@ -0,0 +1,97 @@ +# include "psphotInternal.h" + +bool psphotFullForceReadout(pmConfig *config, const pmFPAview *view, const char *filerule) { + + // measure the total elapsed time in psphotReadout. + psTimerStart ("psphotReadout"); + + // allow objects to be fit with ugly models (central holes, extreme asymmetry, etc) + pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + // set the photcode for this image + if (!psphotAddPhotcode (config, view, filerule)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } + + // optional break-point for processing + char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); + PS_ASSERT_PTR_NON_NULL (breakPt, false); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndVariance (config, view, filerule); + if (!strcasecmp (breakPt, "NOTHING")) { + return psphotReadoutCleanup (config, view, filerule); + } + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + if (!psphotSubtractBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + if (!strcasecmp (breakPt, "BACKMDL")) { + return psphotReadoutCleanup (config, view, filerule); + } + + if (!psphotLoadPSF (config, view, filerule)) { + // this only happens if we had a programming error in psphotLoadPSF + psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); + return psphotReadoutCleanup (config, view, filerule); + } + + // include externally-supplied sources + psphotLoadExtSources (config, view, filerule); + + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + psphotMergeSources (config, view, filerule); + + // generate a psf model for any readouts which need one + psphotFullForcePSF (config, view, filerule); + + // Construct an initial model for each object, set the radius to fitRadius, set circular + // fit mask. NOTE: only applied to sources without guess models + // XXX this currently only generates a PSF model + psphotGuessModels (config, view, filerule); + + // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) + // XXX this will fit non-PSF sources if they are supplied + psphotFitSourcesLinear (config, view, filerule, false, false); + + // measure moments + // XXX does this set the correct gaussian window (how to do multiple windows?) + psphotSourceStats (config, view, filerule, false); + + // measure kron fluxes + psphotKronFlux (config, view, filerule, false); + + // measure petro fluxes + psphotPetroFlux (config, view, filerule); + + // measure radial apertures + psphotRadialApertures (config, view, filerule, 0); + + // measure galaxy shapes + psphotGalaxyShape (config, view, filerule); + + // calculate source magnitudes (psf mag and ap mag) + psphotMagnitudes(config, view, filerule); + + // replace background in residual image + psphotSkyReplace (config, view, filerule); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (config, view, filerule); + + // create the exported-metadata and free local data + return psphotReadoutCleanup (config, view, filerule); +} 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 36146) @@ -0,0 +1,279 @@ +# include "psphotInternal.h" + +bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Fluxes ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // 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->newSources ? detections->newSources : detections->allSources; + psAssert (sources, "missing sources?"); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + // psAssert (psf, "missing psf?"); + + if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) { + + bool status = false; + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping kron fluxes"); + return true; + } + + psTimerStart ("psphot.kron"); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + 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; + } + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } + + // threaded measurement of the source magnitudes + // 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_KRON_ITERATE"); + + 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); + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } +# else + if (!psphotGalaxyShape_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); + return false; + } + 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."); + 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"); + } + psFree(job); + } + } + psFree (cellGroups); + + psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); + return true; +} + +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); + + 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]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + + // check status of this source's moments + if (!source->moments) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + + // XXX where are we storing the supplied kron radius? + if (!isfinite(source->moments->Mrf)) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // XXX this is not right (use same concept as general source fit) + float windowRadius = 6.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); + psAssert (source->pixels, "WTF?"); + + // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr + // do the following for a set of shapes (Ex,Ey) + psphotGalaxyShapeSource (source, maskVal); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + return true; +} + +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) { + + 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; +} + +// 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) { + + 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 + pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize); + + float YY = 0.0; + float YM = 0.0; + float MM = 0.0; + bool usePoisson = false; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[iy][ix] == 0) { + continue; + } + // skip nan value points + if (!isfinite(source->pixels->data.F32[iy][ix])) { + continue; + } + + float fy = source->pixels->data.F32[iy][ix]; + float fm = source->modelFlux->data.F32[iy][ix]; + float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0; + + YY += PS_SQR(fy) * wt; + YM += fm * fy * wt; + MM += PS_SQR(fm) * wt; + } + } + + float Io = YM / MM; + float Chisq = YY - 2 * Io * YM + Io * Io * MM; + + return true; +} Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotImageLoop.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotImageLoop.c (revision 36145) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotImageLoop.c (revision 36146) @@ -132,4 +132,11 @@ } break; + case PSPHOT_FULL_FORCE: + if (!psphotFullForceReadout (config, view, "PSPHOT.INPUT")) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + break; case PSPHOT_MAKE_PSF: if (!psphotMakePSFReadout (config, view, "PSPHOT.INPUT")) { 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 36146) @@ -0,0 +1,259 @@ +# include "psphotInternal.h" + +bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Fluxes ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // 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->newSources ? detections->newSources : detections->allSources; + 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)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) { + + bool status = false; + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping kron fluxes"); + return true; + } + + psTimerStart ("psphot.kron"); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + 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; + } + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } + + // threaded measurement of the source magnitudes + // 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_KRON_ITERATE"); + + 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); + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } +# else + if (!psphotKronFlux_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); + return false; + } + 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."); + 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"); + } + psFree(job); + } + } + psFree (cellGroups); + + psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); + return true; +} + +bool psphotKronFlux_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); + + 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]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + + // check status of this source's moments + if (!source->moments) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + + // XXX where are we storing the supplied kron radius? + if (!isfinite(source->moments->Mrf)) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // is the right?? + float windowRadius = 6.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); + psAssert (source->pixels, "WTF?"); + + // this function populates moments->Mrf,KronFlux,KronFluxErr + psphotKronFluxSource (source, maskVal); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + return true; +} + +// measure just the flux in the given aperture (this is probably a single common function -- can I use one of the pmSourcePhotometry functions?) +bool psphotKronFluxSource (pmSource *source, psImageMaskType maskVal) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + + // the peak position is less accurate but less subject to extreme deviations + float dX = source->moments->Mx - source->peak->xf; + float dY = source->moments->My - source->peak->yf; + float dR = hypot(dX, dY); + float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf; + float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + + // Calculate the Kron magnitude (make this block optional?) + // XXX set the aperture here + float radKron = 2.5*Mrf; + float radKron2 = radKron*radKron; + + int nKronPix = 0; + float Sum = 0.0; + float Var = 0.0; + + // 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; + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radKron) continue; + + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + if (isnan(vPix[row][col])) continue; + + psF32 xDiff = col - xCM; + if (fabs(xDiff) > radKron) continue; + + // radKron is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > radKron2) continue; + + float pDiff = vPix[row][col]; + psF32 wDiff = vWgt[row][col]; + + Sum += pDiff; + Var += wDiff; + nKronPix ++; + } + } + + // return the flux and error to parameters? + source->moments->KronFlux = Sum; + source->moments->KronFluxErr = sqrt(Var); + + return true; +} Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotLoadPSF.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotLoadPSF.c (revision 36145) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotLoadPSF.c (revision 36146) @@ -1,3 +1,20 @@ # include "psphotInternal.h" + +// PSPHOT.PSF.LOAD vs input file -- see note at top +bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule) { + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // Generate the mask and weight images, including the user-defined analysis region of interest + if (!psphotLoadPSFReadout (config, view, filerule, "PSPHOT.PSF.LOAD", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i); + return false; + } + } + return true; +} // NOTE : pmPSF_IO.c functions must load the psf model onto the chip->analysis metadata because @@ -58,19 +75,2 @@ return true; } - -// PSPHOT.PSF.LOAD vs input file -- see note at top -bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule) { - - int num = psphotFileruleCount(config, filerule); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - - // Generate the mask and weight images, including the user-defined analysis region of interest - if (!psphotLoadPSFReadout (config, view, filerule, "PSPHOT.PSF.LOAD", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i); - return false; - } - } - return true; -} Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c (revision 36145) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c (revision 36146) @@ -65,4 +65,7 @@ // Merge the externally supplied sources with the existing sources. Mark them as having mode // PM_SOURCE_MODE_EXTERNAL. + +// XXX this function needs to be updated slightly for psphotFullForce: +// * load the additional parameters to guide the new concepts // XXX This function needs to be updated to loop over set of input files. At the moment, we 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 36146) @@ -0,0 +1,259 @@ +# include "psphotInternal.h" + +bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Petro Fluxes ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // 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->newSources ? detections->newSources : detections->allSources; + 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)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) { + + bool status = false; + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping kron fluxes"); + return true; + } + + psTimerStart ("psphot.kron"); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + 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; + } + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } + + // threaded measurement of the source magnitudes + // 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_KRON_ITERATE"); + + 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); + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } +# else + if (!psphotPetroFlux_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); + return false; + } + 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."); + 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"); + } + psFree(job); + } + } + psFree (cellGroups); + + psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); + return true; +} + +bool psphotPetroFlux_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); + + 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]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + + // check status of this source's moments + if (!source->moments) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + + // XXX where are we storing the supplied kron radius? + if (!isfinite(source->moments->Mrf)) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // is the right?? + float windowRadius = 6.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); + psAssert (source->pixels, "WTF?"); + + // this function populates moments->Mrf,PetroFlux,PetroFluxErr + psphotPetroFluxSource (source, maskVal); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + return true; +} + +// measure just the flux in the given aperture (this is probably a single common function -- can I use one of the pmSourcePhotometry functions?) +bool psphotPetroFluxSource (pmSource *source, psImageMaskType maskVal) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + + // the peak position is less accurate but less subject to extreme deviations + float dX = source->moments->Mx - source->peak->xf; + float dY = source->moments->My - source->peak->yf; + float dR = hypot(dX, dY); + float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf; + float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + + // Calculate the Petro magnitude (make this block optional?) + // XXX set the aperture here + float radPetro = 2.5*Mrf; + float radPetro2 = radPetro*radPetro; + + int nPetroPix = 0; + float Sum = 0.0; + float Var = 0.0; + + // 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; + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radPetro) continue; + + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + if (isnan(vPix[row][col])) continue; + + psF32 xDiff = col - xCM; + if (fabs(xDiff) > radPetro) continue; + + // radPetro is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > radPetro2) continue; + + float pDiff = vPix[row][col]; + psF32 wDiff = vWgt[row][col]; + + Sum += pDiff; + Var += wDiff; + nPetroPix ++; + } + } + + // return the flux and error to parameters? + source->moments->PetroFlux = Sum; + source->moments->PetroFluxErr = sqrt(Var); + + return true; +}