Index: /branches/eam_branch_20080511/ppSim/notes.txt =================================================================== --- /branches/eam_branch_20080511/ppSim/notes.txt (revision 17701) +++ /branches/eam_branch_20080511/ppSim/notes.txt (revision 17702) @@ -1,2 +1,15 @@ + +2008.05.14 + + ppSim needs some work to be a bit more realistic. It currently does + not do the correct thing for generating flats or shutter + corrections. It is only applying them to the background sky model. + It is not applying them to the stellar flux. Also, for adding + sources on top of existing images, it needs to apply the Poisson + stats to the generated stars and then add the pixel values to the + input image. currently, it is placing them directly in the input + image. + + ppSim development work still needed: Index: /branches/eam_branch_20080511/ppSim/src/Makefile.am =================================================================== --- /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17702) @@ -26,4 +26,5 @@ ppSimLoadSpots.c \ ppSimPhotom.c \ + ppSimMergeReadouts.c \ ppSimMosaicChip.c \ ppSimBadPixels.c Index: /branches/eam_branch_20080511/ppSim/src/ppSim.h =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17702) @@ -95,5 +95,5 @@ psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config); -psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) ; +psVector *ppSimMakeBias (bool *status, pmReadout *readout, pmConfig *config, const psRandom *rng) ; bool ppSimMakeDark (pmReadout *readout, pmConfig *config); bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config); @@ -174,3 +174,5 @@ float ppSimGetZeroPoint (psMetadata *recipe, char *filter); +bool ppSimMergeReadouts (pmReadout *outReadout, pmReadout *inReadout); + #endif Index: /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 17702) @@ -40,3 +40,2 @@ return true; } - Index: /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 17702) @@ -4,4 +4,7 @@ bool mdok; + + // skip this step if we did not generate a bias level + if (biasCols == NULL) return true; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe Index: /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 17702) @@ -147,5 +147,5 @@ char *typeStr = ppSimArgToRecipeStr (&status, options, "IMAGE.TYPE", arguments, "-type"); // Requested exposure type - if (!status) { + if (typeStr == NULL) { psError(PS_ERR_IO, false, "An exposure type must be specified using -type"); psFree(arguments); @@ -178,4 +178,5 @@ ppSimArgToRecipeS32(&status, options, "BINNING", arguments, "-bin"); + // XX don't do this: recipe is not yet fully realized (camera options not interpolated) psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe Index: /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17702) @@ -17,12 +17,7 @@ return NULL; } - - simImage = true; - 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; - } } 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) { @@ -30,5 +25,12 @@ return NULL; } - fpa = psMemIncrRefCounter (input->fpa); + } + + // 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; } @@ -98,5 +100,5 @@ // PPSIM.SOURCES carries the constructed, fake sources with their true parameters - // XXX only invoke this code for OBJECT types of images + // XXX only invoke this code for OBJECT types of images? pmFPAfile *simSources = pmFPAfileDefineOutput (config, output->fpa, "PPSIM.SOURCES"); if (!simSources) { @@ -106,6 +108,6 @@ simSources->save = true; - // if we have loaded an input image, we do not need to populate the fpa - if (!simImage) { + // 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... @@ -118,6 +120,4 @@ float zp = ppSimGetZeroPoint (recipe, filter); psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); - - return output; } Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17702) @@ -8,4 +8,6 @@ bool ppSimLoop(pmConfig *config) { + bool status; + PS_ASSERT_PTR_NON_NULL(config, PS_EXIT_PROG_ERROR); @@ -13,4 +15,8 @@ pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, OUTPUT_FILE); // 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 @@ -105,4 +111,9 @@ psVector *biasCols = ppSimMakeBiassec (cell, config); + + pmCell *inputCell = NULL; + if (input) { + inputCell = pmFPAviewThisCell (view, input->fpa); + } for (int i = 0; i < cell->readouts->n; i++) { @@ -125,19 +136,20 @@ } - psVector *biasRows = ppSimMakeBias (readout, config, rng); + psVector *biasRows = ppSimMakeBias (&status, readout, config, rng); + if (!status) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure"); if (type == PPSIM_TYPE_BIAS) goto done; - ppSimMakeDark (readout, config); + if (!ppSimMakeDark (readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure"); if (type == PPSIM_TYPE_DARK) goto done; - ppSimMakeSky (readout, expCorr, type, config); + if (!ppSimMakeSky (readout, expCorr, type, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating sky background"); if (type == PPSIM_TYPE_FLAT) goto done; if (type == PPSIM_TYPE_OBJECT) { - ppSimInsertStars (readout, expCorr, stars, config); + if (!ppSimInsertStars (readout, expCorr, stars, config)) ESCAPE (PS_ERR_UNKNOWN, "problem inserting stars"); } if (type == PPSIM_TYPE_OBJECT) { - ppSimInsertGalaxies (readout, expCorr, galaxies, config); + if (!ppSimInsertGalaxies (readout, expCorr, galaxies, config)) ESCAPE (PS_ERR_UNKNOWN, "problem inserting galaxies"); } @@ -145,10 +157,10 @@ done: - ppSimAddNoise(readout->image, readout->weight, cell, config, rng); - ppSimSaturate(readout, config); - - ppSimBadPixels(readout, config, rng); - - ppSimAddOverscan (readout, config, biasCols, biasRows, rng); + if (!ppSimAddNoise(readout->image, readout->weight, cell, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding noise"); + if (!ppSimSaturate(readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem setting saturation levels"); + + if (!ppSimBadPixels(readout, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding bad pixels"); + + if (!ppSimAddOverscan (readout, config, biasCols, biasRows, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding overscan region"); psFree(biasRows); @@ -156,11 +168,16 @@ readout->parent->data_exists = true; readout->parent->parent->data_exists = true; + + if (inputCell && inputCell->readouts->data[i]) { + if (!ppSimMergeReadouts (readout, inputCell->readouts->data[i])) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image"); + } } psFree(biasCols); - ppSimUpdateConceptsCell (cell, config); + if (!ppSimUpdateConceptsCell (cell, config)) ESCAPE (PS_ERR_UNKNOWN, "problem updating cell concepts"); if (cell->hdu) { - ppSimInitHeader(config, NULL, NULL, cell); + // XXX only do this if there is no INPUT image? + if (!ppSimInitHeader(config, NULL, NULL, cell)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); } @@ -175,9 +192,10 @@ if (chip->hdu) { - ppSimInitHeader(config, NULL, chip, NULL); + // XXX only do this if there is no INPUT image + if (!ppSimInitHeader(config, NULL, chip, NULL)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); } // we perform photometry on the readouts of this chip in the output - ppSimPhotom (config, view); + if (!ppSimPhotom (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem performing photometry"); if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { @@ -185,5 +203,4 @@ psFree(rng); psFree(view); - // return PS_EXIT_SYS_ERROR; return false; } @@ -195,5 +212,6 @@ if (fpa->hdu) { - ppSimInitHeader(config, fpa, NULL, NULL); + // XXX only do this if there is no INPUT image + if (!ppSimInitHeader(config, fpa, NULL, NULL)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); } @@ -202,5 +220,4 @@ psFree(rng); psFree(view); - // return PS_EXIT_SYS_ERROR; return false; } Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c (revision 17702) @@ -1,7 +1,9 @@ # include "ppSim.h" -psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) { +psVector *ppSimMakeBias (bool *status, pmReadout *readout, pmConfig *config, const psRandom *rng) { bool mdok; + + if (status) *status = true; pmCell *cell = readout->parent; @@ -9,6 +11,6 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe - bool bias = psMetadataLookupBool(&status, recipe, "BIAS"); // Generate a Bias? - if (!bias) return true; + bool bias = psMetadataLookupBool(&mdok, recipe, "BIAS"); // Generate a Bias? + if (!bias) return NULL; float biasLevel = psMetadataLookupF32(NULL, recipe, "BIAS.LEVEL"); // Bias level @@ -22,4 +24,5 @@ if (!mdok) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + *status = false; return NULL; } Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c (revision 17702) @@ -11,9 +11,20 @@ psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe - bool dark = psMetadataLookupBool(&status, recipe, "DARK"); // Generate a DARK? + bool dark = psMetadataLookupBool(&mdok, recipe, "DARK"); // Generate a DARK? if (!dark) return true; - float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate - float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + float darkRate = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate + if (isnan(darkRate)) { + psError(PS_ERR_UNKNOWN, false, "missing DARK.RATE\n"); + return false; + } + + float expTime = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time + if (isnan(expTime)) { + psError(PS_ERR_UNKNOWN, false, "missing EXPTIME\n"); + return false; + } + + psTrace("ppSim", 6, "darkRate: %f, expTime: %f\n", darkRate, expTime); for (int y = 0; y < signal->numRows; y++) { Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 17701) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 17702) @@ -1,5 +1,7 @@ # include "ppSim.h" -// XXX add bounds to the inputs? +// this function sets the skyRate to a value for the night sky (SKY.RATE or SKY.MAGS) or for a +// flat-field image (FLAT.RATE). Include a shutter correction and a scattered light source + bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config) { @@ -16,16 +18,20 @@ bool sky = psMetadataLookupBool(&status, recipe, "SKY"); // Generate a SKY flux? - if (!sky) return true; - float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time - float flatSigma = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient - float flatRate = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate - float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time - float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate + float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time + + float flatSigma = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient + float flatRate = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate + float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time + float scatterScale = psMetadataLookupF32(&status, recipe, "SCATTER.SCALE"); // scattered light term + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { - float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); - float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status); - float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status); + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + if (type == PPSIM_TYPE_FLAT) { + skyRate = flatRate; } @@ -44,8 +50,10 @@ // Size of FPA psRegion *bounds = ppSimFPABounds (fpa); + int dXfpa = bounds->x1 - bounds->x0; + int dYfpa = bounds->y1 - bounds->y0; // Correct chip offsets so that boresight is in the middle of the FPA - x0Chip -= 0.5 * (bounds->x1 - bounds->x0); - y0Chip -= 0.5 * (bounds->y1 - bounds->y0); + x0Chip -= 0.5 * dXfpa; + y0Chip -= 0.5 * dYfpa; for (int y = 0; y < signal->numRows; y++) { @@ -62,20 +70,25 @@ // Gaussian flat-field over the FPA - float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) / - flatSigma / sqrtf(2.0 * M_PI); + float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) / flatSigma / sqrtf(2.0 * M_PI); - if (type == PPSIM_TYPE_FLAT) { - float flatFlux = flatRate * flatValue * realExpTime; // Flux from flat-field - signal->data.F32[y][x] += flatFlux; - variance->data.F32[y][x] += flatFlux; - continue; - } + float scatterRate = 0.0; - expCorr->data.F32[y][x] = realExpTime / expTime; + if (sky) { + // add a scattered light term to the flat-field images (no + if (type == PPSIM_TYPE_FLAT) { + float xF = 2.0*(xFPA / dXfpa) - 1.0; + scatterRate = scatterScale * PS_SQR(xF); + } - // Sky background - float skyFlux = skyRate * flatValue * realExpTime; // Flux from sky - signal->data.F32[y][x] += skyFlux; - variance->data.F32[y][x] += skyFlux; + // Sky background + float skyFlux = (skyRate * flatValue + scatterRate) * realExpTime; // Flux from sky + signal->data.F32[y][x] += skyFlux; + variance->data.F32[y][x] += skyFlux; + } + + // used later to modify the star and galaxy photometry + if (expCorr) { + expCorr->data.F32[y][x] = realExpTime / expTime; + } // TO DO: Add fringes