Index: trunk/ppSim/src/ppSimCreate.c =================================================================== --- trunk/ppSim/src/ppSimCreate.c (revision 17915) +++ trunk/ppSim/src/ppSimCreate.c (revision 18011) @@ -1,127 +1,226 @@ # include "ppSim.h" - -// XXX this function forces us to define the camera (on the command line) and format (via the -// PPSIM.OUTPUT entry). In this case, we need to set config->format,formatName based on these -// values. This will be a problem when we want to load an input image (in order to add fake -// stars). We will need to add some logic in ppSimArguments to distinguish the cases of 1) -// input image and 2) specified camera pmFPAfile *ppSimCreate(pmConfig *config) { + PS_ASSERT_PTR_NON_NULL(config, 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"); + 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; - } - - simImage = true; - fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the observation - if (!fpa) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); - return NULL; - } - + // 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 { - simImage = false; - if (input->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_IO, true, "PPIMAGE.INPUT is not of type IMAGE"); - return NULL; - } - fpa = input->fpa; - } - - // define the output image file - pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE); - if (!file) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s. " - "Did you forget to specify the format?", OUTPUT_FILE); - return NULL; - } - if (file->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "%s type is not IMAGE", OUTPUT_FILE); - psFree(fpa); - psFree(file); - return NULL; - } - file->save = true; - - config->format = psMemIncrRefCounter (file->format); - config->formatName = psStringCopy (file->formatName); - - // have we supplied a psf model? - if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { - bool status = false; - - // tie the psf file to the chipMosaic - pmFPAfileBindFromArgs(&status, file, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); - psFree(fpa); - psFree(file); - return NULL; - } - } - - // XXX only invoke this code for OBJECT types of images + // If an image is supplied, we still generate a fake image and merge them together downstream + // (otherwise, we get the variance wrong). + 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, config->cameraName); // 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 we have been supplied an input image, but no ra & dec, use the input image values + if (input) { + float ra = psMetadataLookupF32(&status, recipe, "RA"); + if (isnan(ra)) { + ra = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.RA"); + psMetadataAddF32(recipe, PS_LIST_TAIL, "RA", PS_META_REPLACE, "ra (radians)", ra); + } + float dec = psMetadataLookupF32(&status, recipe, "DEC"); + if (isnan(dec)) { + dec = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.DEC"); + psMetadataAddF32(recipe, PS_LIST_TAIL, "DEC", PS_META_REPLACE, "dec (radians)", dec); + } + } + } + + 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) { + + // XXX at the moment, we can perform photometry on the fake sources and the forced + // photometry positions, but not one or the other. Also, we only support photometry in + // the context of an image previously analysed by psphot. Add support for: + // * photometry of a pure fake image (requires peak detection and psf creation) + // * photometry of forced source positions without fake image (this might work, it just + // requires not generating any fake features). + + if (!input) { + psError(PS_ERR_UNKNOWN, false, "input image not found; currently required for photometry"); + return NULL; + } + + // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) + pmFPAfile *fakeImage = pmFPAfileDefineChipMosaic(config, output->fpa, "PPSIM.FAKE.CHIP"); + if (!fakeImage) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FAKE.CHIP")); + psFree(fpa); + return NULL; + } + if (fakeImage->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.FAKE.CHIP is not of type IMAGE"); + psFree(fpa); + return NULL; + } + + // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) + pmFPAfile *forceImage = NULL; + if (input) { + forceImage = pmFPAfileDefineChipMosaic(config, input->fpa, "PPSIM.FORCE.CHIP"); + if (!forceImage) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FORCE.CHIP")); + psFree(fpa); + return NULL; + } + if (forceImage->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.FORCE.CHIP is not of type IMAGE"); + psFree(fpa); + return NULL; + } + } + + // define associated psphot input/output files + if (!ppSimPhotomFiles (config, fakeImage, forceImage)) { + 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 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; + } + } + } + // PPSIM.SOURCES carries the constructed, fake sources with their true parameters - pmFPAfile *simSources = pmFPAfileDefineOutput (config, file->fpa, "PPSIM.SOURCES"); + // 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; + 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 do not need to populate the fpa - if (!simImage) { - return file; - } - - pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level at which PHU goes + // 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, file->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - psFree(file); - return NULL; - } + 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, file->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - psFree(file); - return NULL; - } - } - - pmCell *cell; // Cell from chip - while ((cell = pmFPAviewNextCell(view, fpa, 1))) { - if (phuLevel == PM_FPA_LEVEL_CELL) { - if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - psFree(file); - return NULL; - } - } - } + 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; + } + } + } } @@ -129,4 +228,4 @@ psFree(view); - return file; + return output; }