Index: /branches/eam_branch_20080511/ppSim/notes.txt =================================================================== --- /branches/eam_branch_20080511/ppSim/notes.txt (revision 17817) +++ /branches/eam_branch_20080511/ppSim/notes.txt (revision 17818) @@ -1,2 +1,41 @@ + + +2008.05.17 + + pmFPAfiles used in ppSim: + + PPSIM.INPUT : an input image loaded by ppSim; the simulated features + may be added to this image. In fact, the features are + generated in the output image, and the two images are + added (as images) together. + + PPSIM.OUTPUT : the output image. this file is used for the main + ppSimLoop function + + PPSIM.CHIP : chip-mosaiced version of the output image. only + generated if photometry is requested. Is this + actually used as an output image? + + PPSIM.REAL.SOURCES : real sources loaded by PPSIM. these sources + are subtracted from the image before testing + for the detectability of the fake sources and + before measuring the flux of fake or forced + soures. They should be the sources measured in + a previous psphot run, though any input list + could be used. Results are very ill-defined if + the sources do not correspond to actual image + sources! + + PSPHOT.PSF.LOAD : input PSF model used to measure the input sources + + PPSIM.SOURCES : output fake sources (injected fake sources with real values) + PPSIM.FAKE.SOURCES : output fake photometry (measured photometry for all fake sources) + PPSIM.FORCE.SOURCES : output force photometry + + ** used within psphot functions: + PSPHOT.BACKMDL : model background generated by psphotModelBackground + PSPHOT.BACKMDL.STDEV : model background error generated by psphotModelBackground + PSPHOT.BACKGND : full-scale model background generated by psphotSubtractBackground + PSPHOT.BACKSUB : background-subtracted image 2008.05.15 Index: /branches/eam_branch_20080511/ppSim/src/Makefile.am =================================================================== --- /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17818) @@ -26,5 +26,10 @@ ppSimLoadSpots.c \ ppSimPhotom.c \ + ppSimPhotomReadout.c \ + ppSimPhotomFiles.c \ + ppSimLoadForceSources.c \ ppSimMergeReadouts.c \ + ppSimDetections.c \ + ppSimMergeSources.c \ ppSimMosaicChip.c \ ppSimRandomGaussian.c \ Index: /branches/eam_branch_20080511/ppSim/src/ppSim.h =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17818) @@ -16,5 +16,5 @@ #define PPSIM_RECIPE "PPSIM" -#define OUTPUT_FILE "PPSIM.OUTPUT" +// #define OUTPUT_FILE "PPSIM.OUTPUT" // Compare a value with minimum and maximum values, replacing where required. @@ -174,5 +174,5 @@ float ppSimGetZeroPoint (psMetadata *recipe, char *filter); -bool ppSimMergeReadouts (pmReadout *outReadout, pmReadout *inReadout); +bool ppSimMergeReadouts (pmConfig *config, pmFPAview *view); double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma); @@ -180,3 +180,9 @@ void ppSimRandomGaussianFree(); +bool ppSimPhotomFiles (pmConfig *config, pmFPAfile *input); +bool ppSimPhotomReadout(pmConfig *config, const pmFPAview *view); +psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view); +bool ppSimDetections (psImage *significance, psMetadata *recipe, psArray *sources); +psArray *ppSimMergeSources (psArray *in1, psArray *in2); + #endif Index: /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17818) @@ -3,204 +3,188 @@ pmFPAfile *ppSimCreate(pmConfig *config) { - bool status; - bool simImage = false; - PS_ASSERT_PTR_NON_NULL(config, NULL); - pmFPA *fpa = NULL; + bool status; + bool simImage = false; + PS_ASSERT_PTR_NON_NULL(config, NULL); + pmFPA *fpa = NULL; - // the input image defines the camera. if it is not supplied, the user must have - // supplied a camera and other metadata on the command line - pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPIMAGE.INPUT", "INPUT"); - if (!input) { - // if we have not specified the camera already, we need to interpolate the recipes associated with this camera, and read other command-line recipes - if (!pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) { - psError(PS_ERR_IO, false, "Error merging recipes from camera config for %s", config->cameraName); - return NULL; - } - } else { - // If an image is supplied, we still generate a fake image and merge them together downstream - // (otherwise, we get the variance wrong). - simImage = false; - if (input->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_IO, true, "PPIMAGE.INPUT is not of type IMAGE"); - return NULL; - } + // the input image defines the camera. if it is not supplied, the user must have + // supplied a camera and other metadata on the command line + pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPSIM.INPUT", "INPUT"); + if (!input) { + // if we have not specified the camera already, we need to interpolate the recipes associated with this camera, and read other command-line recipes + if (!pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) { + psError(PS_ERR_IO, false, "Error merging recipes from camera config for %s", config->cameraName); + return NULL; } + } else { + // If an image is supplied, we still generate a fake image and merge them together downstream + // (otherwise, we get the variance wrong). + simImage = false; + if (input->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.INPUT is not of type IMAGE"); + return NULL; + } + } - // generate the fpa structure used by the output camera (determined from INPUT or specified) - assert (config->camera); - fpa = pmFPAConstruct(config->camera); // FPA to contain the observation - if (!fpa) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + // generate the fpa structure used by the output camera (determined from INPUT or specified) + assert (config->camera); + fpa = pmFPAConstruct(config->camera); // FPA to contain the observation + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + return NULL; + } + + // define the output image file -- this is the basis for the ppSimLoop + pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSIM.OUTPUT"); + if (!output) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from PPSIM.OUTPUT. Did you forget to specify the format?"); + return NULL; + } + if (output->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "PPSIM.OUTPUT type is not IMAGE"); + psFree(fpa); + return NULL; + } + // XXX we should not require the output image to be written + output->save = true; + + config->format = psMemIncrRefCounter (output->format); + config->formatName = psStringCopy (output->formatName); + + // the recipe is now fully realized for the desired camera + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate + + if (type == PPSIM_TYPE_OBJECT) { + // adjust the seeing by the scale + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); + if (isnan(seeing)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "seeing is not defined"); + psFree(fpa); + return NULL; + } + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); + if (isnan(scale)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "pixel scale is not defined"); + psFree(fpa); + return NULL; + } + psMetadataAddF32(recipe, PS_LIST_TAIL, "SEEING", PS_META_REPLACE, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); + } + + if ((type == PPSIM_TYPE_OBJECT) || (type == PPSIM_TYPE_FLAT)) { + // determine the zeropoint from the filter + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); + if (isnan(zp)) { + char *filter = psMetadataLookupStr(&status, recipe, "FILTER"); + float zp = ppSimGetZeroPoint (recipe, filter); + psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + } + } + + // For photometry, we operate on the chip-mosaicked image. we create a copy of the mosaicked + // image for psphot so we can write out a clean image + bool doPhotom = psMetadataLookupBool(&status, recipe, "PHOTOM"); // Density of fakes + if (doPhotom) { + + // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) + pmFPAfile *chipImage = pmFPAfileDefineChipMosaic(config, output->fpa, "PPSIM.CHIP"); + if (!chipImage) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.CHIP")); + psFree(fpa); + return NULL; + } + if (chipImage->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.CHIP is not of type IMAGE"); + psFree(fpa); return NULL; } - // define the output image file -- this is the basis for the ppSimLoop - pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE); - if (!output) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s. " - "Did you forget to specify the format?", OUTPUT_FILE); - return NULL; + // define associated psphot input/output files + if (!ppSimPhotomFiles (config, chipImage)) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files for psphot"); + psFree(fpa); + return NULL; } - if (output->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "%s type is not IMAGE", OUTPUT_FILE); + } else { + // have we supplied a psf model? this happens in ppSimPhotomFiles if we request a photometry + // analysis. however, even if we do not, a psf model may be used to generate the fake + // sources. + if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { + // tie the psf file to the chipMosaic + pmFPAfileBindFromArgs(&status, output, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); psFree(fpa); - return NULL; + return NULL; + } } - // XXX we should not require the output image to be written - output->save = true; + } - config->format = psMemIncrRefCounter (output->format); - config->formatName = psStringCopy (output->formatName); + // PPSIM.SOURCES carries the constructed, fake sources with their true parameters + // XXX only invoke this code for OBJECT types of images? + pmFPAfile *simSources = pmFPAfileDefineOutput (config, output->fpa, "PPSIM.SOURCES"); + if (!simSources) { + psError(PS_ERR_UNKNOWN, false, "Cannot find a rule for PPSIM.SOURCES"); + return false; + } + simSources->save = true; - // the recipe is now fully realized for the desired camera + // if we have loaded an input image, we derive certain values from the image, if possible + if (input) { + // we need to extract certain metadata from the image and populate the recipe. + // or else we need to set the fpa concepts based on the recipe options... + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate - ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate + ppSimArgToRecipeF32(&status, recipe, "EXPTIME", input->fpa->concepts, "FPA.EXPOSURE"); + char *filter = ppSimArgToRecipeStr(&status, recipe, "FILTER", input->fpa->concepts, "FPA.FILTER"); - if (type == PPSIM_TYPE_OBJECT) { - // adjust the seeing by the scale - float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); - if (isnan(seeing)) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "seeing is not defined"); - psFree(fpa); - return NULL; - } - float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); - if (isnan(scale)) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "pixel scale is not defined"); - psFree(fpa); - return NULL; - } - psMetadataAddF32(recipe, PS_LIST_TAIL, "SEEING", PS_META_REPLACE, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); + float zp = ppSimGetZeroPoint (recipe, filter); + psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + } + + pmFPALevel phuLevel = pmFPAPHULevel(output->format); // Level at which PHU goes + + pmFPAview *view = pmFPAviewAlloc(0);// View for current level + + if (phuLevel == PM_FPA_LEVEL_FPA) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } + } + + pmChip *chip; // Chip from FPA + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + if (phuLevel == PM_FPA_LEVEL_CHIP) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } } - if ((type == PPSIM_TYPE_OBJECT) || (type == PPSIM_TYPE_FLAT)) { - // determine the zeropoint from the filter - float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); - if (isnan(zp)) { - char *filter = psMetadataLookupStr(&status, recipe, "FILTER"); - float zp = ppSimGetZeroPoint (recipe, filter); - psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + pmCell *cell; // Cell from chip + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + if (phuLevel == PM_FPA_LEVEL_CELL) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; } + } } + } - // For photometry, we operate on the chip-mosaicked image. we create a copy of the mosaicked - // image for psphot so we can write out a clean image - bool doPhotom = psMetadataLookupBool(&status, recipe, "PHOTOM"); // Density of fakes - if (doPhotom) { + psFree(fpa); + psFree(view); - // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) - pmFPAfile *chipImage = pmFPAfileDefineChipMosaic(config, output->fpa, "PPSIM.CHIP"); - if (!chipImage) { - psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.CHIP")); - psFree(fpa); - return NULL; - } - if (chipImage->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_IO, true, "PPSIM.CHIP is not of type IMAGE"); - psFree(fpa); - return NULL; - } - - // *** Define two input files: - // PPSIM.REAL.SOURCES - - // the input image(s) are required arguments; they define the camera - pmFPAfile *realSources = pmFPAfileDefineFromArgs (&status, config, "PPSIM.REAL.SOURCES", "REAL.SOURCES"); - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to build FPA from PPSIM.REAL.SOURCES"); - return false; - } - - // this file is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by - // psphotDefineFiles - pmFPAfile *psphotInput = pmFPAfileDefineFromFile (config, chipImage, 1, 1, "PSPHOT.INPUT"); - PS_ASSERT (psphotInput, false); - - // define associated psphot input/output files - if (!psphotDefineFiles (config, psphotInput)) { - psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files for psphot"); - psFree(fpa); - return NULL; - } - } else { - // have we supplied a psf model? this happens in psphotDefineFiles if we request a photometry - // analysis. however, even if we do not, a psf model may be used to generate the fake - // sources. - if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { - // tie the psf file to the chipMosaic - pmFPAfileBindFromArgs(&status, output, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); - psFree(fpa); - return NULL; - } - } - } - - // PPSIM.SOURCES carries the constructed, fake sources with their true parameters - // XXX only invoke this code for OBJECT types of images? - pmFPAfile *simSources = pmFPAfileDefineOutput (config, output->fpa, "PPSIM.SOURCES"); - if (!simSources) { - psError(PS_ERR_UNKNOWN, false, "Cannot find a rule for PPSIM.SOURCES"); - return false; - } - simSources->save = true; - - // if we have loaded an input image, we derive certain values from the image, if possible - if (input) { - // we need to extract certain metadata from the image and populate the recipe. - // or else we need to set the fpa concepts based on the recipe options... - - psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - - ppSimArgToRecipeF32(&status, recipe, "EXPTIME", input->fpa->concepts, "FPA.EXPOSURE"); - char *filter = ppSimArgToRecipeStr(&status, recipe, "FILTER", input->fpa->concepts, "FPA.FILTER"); - - float zp = ppSimGetZeroPoint (recipe, filter); - psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); - } - - pmFPALevel phuLevel = pmFPAPHULevel(output->format); // Level at which PHU goes - - pmFPAview *view = pmFPAviewAlloc(0);// View for current level - - if (phuLevel == PM_FPA_LEVEL_FPA) { - if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - return NULL; - } - } - - pmChip *chip; // Chip from FPA - while ((chip = pmFPAviewNextChip(view, fpa, 1))) { - if (phuLevel == PM_FPA_LEVEL_CHIP) { - if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - return NULL; - } - } - - pmCell *cell; // Cell from chip - while ((cell = pmFPAviewNextCell(view, fpa, 1))) { - if (phuLevel == PM_FPA_LEVEL_CELL) { - if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - return NULL; - } - } - } - } - - psFree(fpa); - psFree(view); - - return output; + return output; } Index: /branches/eam_branch_20080511/ppSim/src/ppSimDetectionLimits.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimDetectionLimits.c (revision 17818) +++ /branches/eam_branch_20080511/ppSim/src/ppSimDetectionLimits.c (revision 17818) @@ -0,0 +1,47 @@ +# include "ppSim.h" + +// compare injected sources (PPSIM.SOURCES) and measured fake sources (PPSIM.FAKE.SOURCES) +// to determine detection limits and recovered magnitude errors +bool ppSimDetectionLimits (pmConfig *config, pmFPAview *view) { + + // 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; + } + + // XXX do we need to ask if the analysis was performed to 1st or 2nd detection limit? + float NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + psArray *injectedSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.SOURCES"); + psArray *measuredSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.FAKE.SOURCES"); + psAssert (injectedSources->n == measuredSources->n, "mis-match between injected and measured sources"); + + psVector *mag = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *dmag = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *SN = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *detect = psVectorAlloc (injectedSources->n, PS_TYPE_BOOL); + + // first measure, for each source, if it was detected (SN > limit) + // and the magnitude offset (dM = M_inject - M_measure) + for (int i = 0; i < injectedSources->n; i++) { + pmSource *injectSource = injectedSources->data[i]; + pmSource *measureSource = measuredSources->data[i]; + + SN->data.F32[i] = measureSource->peak->SN; + dmag->data.F32[i] = injectSource->psfMag - measureSource->psfMag; + mag->data.F32[i] = injectSource->psfMag; + detect->data.Bool[i] = (measureSource->peak->SN >= NSIGMA_PEAK); + } + + // generate a histogram consisting of the instrumental magnitude bin, the + // XXX what is resolution? (user parameter?) + + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimDetections.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimDetections.c (revision 17818) +++ /branches/eam_branch_20080511/ppSim/src/ppSimDetections.c (revision 17818) @@ -0,0 +1,29 @@ +# include "ppSim.h" + +bool ppSimDetections (psImage *significance, psMetadata *recipe, psArray *sources) { + + bool status; + + // for each source, measure the significance of the peak at the given coordinate + // where does this get stored? + + // XXX need to get the effective Area from the PSF sigma + float SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "SIGMA_SMOOTH"); + psAssert (status, "SIGMA_SMOOTH missing: call psphotSignificanceImage first"); + float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH); + + int row0 = significance->row0; + int col0 = significance->col0; + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + pmPeak *peak = source->peak; + psAssert (peak, "peak is not defined for the source"); + + peak->value = significance->data.F32[peak->y-row0][peak->x-col0]; + peak->SN = sqrt(peak->value*effArea); + + } + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 17818) @@ -15,7 +15,4 @@ if (!stars->n) { return true; } - - // XXX is this needed? - // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources pmCell *cell = readout->parent; @@ -136,5 +133,6 @@ // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT" - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); + // XXX PPSIM.SOURCES are the fake sources with their *input* values + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PPSIM.SOURCES", PS_DATA_ARRAY, "fake sources", sources); psFree(sources); Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c (revision 17818) @@ -1,2 +1,8 @@ +# include "ppSim.h" + +psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view) { + +# if (0) + // de-activate all files except PSASTRO.REFSTARS pmFPAfileActivate (config->files, false, NULL); @@ -77,2 +83,6 @@ pmFPAfileActivate (config->files, true, NULL); pmFPAfileActivate (config->files, false, "PSASTRO.OUT.REFSTARS"); +# endif + + return NULL; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17818) @@ -13,10 +13,6 @@ // in this program, we are looping over the output image, rather than the input as in ppImage - pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, OUTPUT_FILE); // Output file + pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); // Output file assert(file); - - // if we have an input image, we need to add the synthetic data on top of it below - // XXX we may potentially use this input file to skip missing elements - pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPIMAGE.INPUT"); pmFPA *fpa = file->fpa; // FPA for file @@ -111,9 +107,4 @@ psVector *biasCols = ppSimMakeBiassec (cell, config); - - pmCell *inputCell = NULL; - if (input) { - inputCell = pmFPAviewThisCell (view, input->fpa); - } for (int i = 0; i < cell->readouts->n; i++) { @@ -169,14 +160,6 @@ readout->parent->parent->data_exists = true; - if (inputCell && inputCell->readouts->data[i]) { - pmReadout *inReadout = inputCell->readouts->data[i]; - if (!inReadout->weight) { - if (!pmReadoutGenerateWeight(inReadout, true)) { - psError (PS_ERR_UNKNOWN, false, "trouble creating weight"); - return false; - } - } - if (!ppSimMergeReadouts (readout, inReadout)) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image"); - } + // if there is an input image, merge it with the simulated image + if (!ppSimMergeReadouts (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image"); } psFree(biasCols); Index: /branches/eam_branch_20080511/ppSim/src/ppSimMergeReadouts.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMergeReadouts.c (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMergeReadouts.c (revision 17818) @@ -2,7 +2,30 @@ // XXX add bounds to the inputs? -bool ppSimMergeReadouts (pmReadout *outReadout, pmReadout *inReadout) { +bool ppSimMergeReadouts (pmConfig *config, pmFPAview *view) { // bool mdok; + + // if we have an input image, we need to add the synthetic data on top of it below + // XXX we may potentially use this input file to skip missing elements + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSIM.INPUT"); + if (!input) return true; + + pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa); + if (!inReadout) return true; + + if (!inReadout->weight) { + if (!pmReadoutGenerateWeight(inReadout, true)) { + psError (PS_ERR_UNKNOWN, false, "trouble creating weight"); + return false; + } + } + + // output must exist or we made a programming error + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); // Output file + assert(output); + + // XXX require outReadout? + pmReadout *outReadout = pmFPAviewThisReadout (view, output->fpa); + if (!outReadout) return true; psImage *inSignal = inReadout->image; @@ -12,5 +35,4 @@ psImage *outVariance = outReadout->weight; - // psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe assert (inSignal->numRows == outSignal->numRows); assert (inSignal->numCols == outSignal->numCols); Index: /branches/eam_branch_20080511/ppSim/src/ppSimMergeSources.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMergeSources.c (revision 17818) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMergeSources.c (revision 17818) @@ -0,0 +1,17 @@ +# include "ppSim.h" + +psArray *ppSimMergeSources (psArray *in1, psArray *in2) { + + psArray *out = psArrayAlloc (in1->n + in2->n); + + for (int i = 0; i < in1->n; i++) { + out->data[i] = psMemIncrRefCounter (in1->data[i]); + } + + int nOff = in1->n; + for (int i = nOff; i < out->n; i++) { + out->data[i] = psMemIncrRefCounter (in2->data[i-nOff]); + } + + return out; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c (revision 17818) @@ -18,13 +18,13 @@ ppSimMosaicChip(config, blankMask, view, "PPSIM.CHIP", "PPSIM.OUTPUT"); - // find or define a pmFPAfile PSPHOT.INPUT - pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + // XXX have ppSimMosaicChip return the pointer? + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PPSIM.CHIP"); if (!status) { - psError(PSPHOT_ERR_CONFIG, false, "PSPHOT.INPUT I/O file is not defined"); + psError(PSPHOT_ERR_CONFIG, false, "PPSIM.CHIP I/O file is not defined"); return false; } // we make a new copy of the output chip to keep psphot from modifying the output image - // XXX is this needed for ppSim? (yes, unless we do not do photometry and make a simulated image) + // XXX is this needed for ppSim? only if we want to write out the image with and without subtracted sources pmChip *oldChip = pmFPAviewThisChip (view, input->src); pmChip *newChip = pmFPAviewThisChip (view, input->fpa); @@ -45,5 +45,4 @@ return false; } - ppSimDetectionLimits (config, view); } Index: /branches/eam_branch_20080511/ppSim/src/ppSimPhotomFiles.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimPhotomFiles.c (revision 17818) +++ /branches/eam_branch_20080511/ppSim/src/ppSimPhotomFiles.c (revision 17818) @@ -0,0 +1,74 @@ +# include "ppSim.h" + +// define the needed / desired I/O files +bool ppSimPhotomFiles (pmConfig *config, pmFPAfile *input) { + + bool status = false; + + // select recipe options supplied on command line + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + + // optionally save the background model (small FITS image) + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) { + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.BACKMDL"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKMDL"); + return false; + } + output->save = true; + } + // optionally save the background model's standard deviation (small FITS image) + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL.STDEV")) { + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.BACKMDL.STDEV"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKMDL.STDEV"); + return false; + } + output->save = true; + } + // optionally save the full background image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKGND")) { + pmFPAfile *output = pmFPAfileDefineFromFile (config, input, 1, 1, "PSPHOT.BACKGND"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKGND"); + return false; + } + output->save = true; + } + // optionally save the background-subtracted image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKSUB")) { + pmFPAfile *output = pmFPAfileDefineFromFile (config, input, 1, 1, "PSPHOT.BACKSUB"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT."); + return false; + } + output->save = true; + } + + // the input sources are required for sensible measurements of the fake or force sources + pmFPAfileDefineFromArgs (&status, config, "PPSIM.REAL.SOURCES", "INPUT.SOURCES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to build FPA from PPSIM.REAL.SOURCES"); + return false; + } + + // the fake sources are carried on the input->fpa structures + pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, input, "PPSIM.FAKE.SOURCES"); + if (!outsources) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.OUTPUT"); + return false; + } + outsources->save = true; + + pmFPAfileBindFromArgs(&status, input, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!status) { + psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.PSF.LOAD"); + return status; + } + + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c (revision 17817) +++ /branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c (revision 17818) @@ -1,5 +1,5 @@ -# include "psphotInternal.h" +# include "ppSim.h" -bool psphotReadout(pmConfig *config, const pmFPAview *view) { +bool ppSimPhotomReadout(pmConfig *config, const pmFPAview *view) { psTimerStart ("psphotReadout"); @@ -18,7 +18,7 @@ } - // find the currently selected readout - // XXX keep this name or go with something like PPSIM.PHOT.INPUT? - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.CHIP"); PS_ASSERT_PTR_NON_NULL (readout, false); @@ -27,4 +27,5 @@ // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + // this function uses PSPHOT.PSF.LOAD as the pmFPAfile pmPSF *psf = psphotLoadPSF (config, view, recipe); assert (psf); @@ -33,15 +34,22 @@ // PPSIM.REAL.SOURCES carries the pmSource objects (from psphot analysis or loaded externally) psArray *realSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.REAL.SOURCES"); - psArray *fakeSources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); - forceSources = ppSimLoadForceSources (config, view); + psArray *forceSources = ppSimLoadForceSources (config, view); - // replace all sources - psphotRemoveSources (realSources, recipe); + // XXX make a copy of the fake-source parameters so these can be saved independently + // from the measured fake-source parameters + psArray *injectedSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.SOURCES"); + psArray *fakeSources = psArrayAlloc (injectedSources->n); + for (int i = 0; i < injectedSources->n; i++) { + fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]); + } + + // remove all sources + psphotRemoveAllSources (realSources, recipe); // generate a background model (median, smoothed image) - if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) { + if (!psphotModelBackground (config, view, "PPSIM.CHIP")) { return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); } - if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { + if (!psphotSubtractBackground (config, view, "PPSIM.CHIP")) { return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); } @@ -51,4 +59,5 @@ // XXX fake sources should measure peak->x,y, force sources should not + psMaskType maskVal = 0xff; psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); ppSimDetections (significance, recipe, fakeSources); @@ -60,5 +69,5 @@ // replace all sources - psphotReplaceAll (realSources, recipe); + psphotReplaceAllSources (realSources, recipe); // construct an initial model for each object @@ -67,14 +76,16 @@ psphotGuessModels (readout, forceSources, recipe, psf); + psArray *sources = NULL; + // linear fit to real + fake sources - psArray *sources = ppSimMergeSources (realSources, fakeSources); + sources = ppSimMergeSources (realSources, fakeSources); psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract - psphotReplaceAll (sources, recipe); + psphotReplaceAllSources (sources, recipe); psFree (sources); // only frees the merged references // linear fit to real + forced sources - psArray *sources = ppSimMergeSources (realSources, forceSources); + sources = ppSimMergeSources (realSources, forceSources); psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract - psphotReplaceAll (sources, recipe); + psphotReplaceAllSources (sources, recipe); psFree (sources); // only frees the merged references @@ -85,15 +96,26 @@ if (!psphotApResid (readout, forceSources, recipe, psf)) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, readout, recipe, NULL, psf, NULL); } - // calculate source magnitudes + // calculate source magnitudes (for which set??) pmReadout *background = psphotSelectBackground (config, view, false); - psphotMagnitudes(sources, recipe, psf, background); + psphotMagnitudes(fakeSources, recipe, psf, background); + psphotMagnitudes(forceSources, recipe, psf, background); // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); + psphotSourceFreePixels (realSources); + psphotSourceFreePixels (fakeSources); + psphotSourceFreePixels (forceSources); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? + // (or don't supply the sources, and do this with a different function) + psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); + + // add forceSources to the readout (real and fake already there) + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PPSIM.FAKE.SOURCES", PS_DATA_ARRAY, "fake photometry ", fakeSources); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PPSIM.FORCE.SOURCES", PS_DATA_ARRAY, "forced photometry ", forceSources); + psFree (forceSources); + return true; }