Index: /branches/eam_branch_20080511/ppSim/notes.txt =================================================================== --- /branches/eam_branch_20080511/ppSim/notes.txt (revision 17708) +++ /branches/eam_branch_20080511/ppSim/notes.txt (revision 17709) @@ -1,13 +1,42 @@ -2008.05.14 +2008.05.15 - 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. + For fake and force photometry, we need to load the known existing + sources, then perform a complete linear fit solution to the complete + set of both real (known) sources plus the fake and/or forced + photometry positions (as PSFs). The fake and forced photometry + positions need to be measured independently. This provides the real + photometry for these sources ** in the presence of the other real + sources **. In order to test the detectability of the fake sources + (or the forced sources, for that matter), we need to generate the + smoothed detection image, and determine the peak value at the + positions of the sources. + + ppSimPhotomReadout outline: + + * we have three source lists (loaded before the function is called): + * realSources : these are loaded from a CMF-style file (or equiv) + * fakeSources : these are defined internally, and need a link from + the input fake source to the measured fake source + * forceSources : these are loaded from a DVO database via getstar + + * we require a supplied PSF + + * need to subtract the background (before or after subtracting the + sources?) + + * subtract the real sources + * model the background + * build the detection image + * measure the detectability of the fake and force sources + * replace the real sources + + * merge real + fake sources + * linear fit to merged source list + * replace the sources (make 'no-subtract' option?) + + * merge real + forced sources + * linear fit to merged source list + * replace the sources (make 'no-subtract' option?) Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c (revision 17709) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c (revision 17709) @@ -0,0 +1,78 @@ + // de-activate all files except PSASTRO.REFSTARS + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSASTRO.OUT.REFSTARS"); + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + // XXX there can only be one readout per chip in astrometry, right? + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, fpa); + + int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1"); + int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2"); + + float minX = -fieldPadding*Nx; + float maxX = (1+fieldPadding)*Nx; + float minY = -fieldPadding*Ny; + float maxY = (1+fieldPadding)*Ny; + + // the refstars is a subset within range of this chip + psArray *refstars = psArrayAllocEmpty (100); + + // select the reference objects within range of this readout + // project the reference objects to this chip + for (int i = 0; i < refs->n; i++) { + pmAstromObj *ref = pmAstromObjCopy(refs->data[index->data.S32[i]]); + + psProject (ref->TP, ref->sky, fpa->toSky); + psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP); + psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); + + // limit the X,Y range of the refs to the selected chip + if (ref->chip->x < minX) goto skip; + if (ref->chip->x > maxX) goto skip; + if (ref->chip->y < minY) goto skip; + if (ref->chip->y > maxY) goto skip; + + psArrayAdd (refstars, 100, ref); + skip: + psFree (ref); + + if (nMax && (refstars->n >= nMax)) break; + } + psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n); + + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars); + psFree (refstars); + + if (matchLumFunc) { + // in this case, no PSASTRO.REFSTARS is added to readout->analysis + if (!psastroRefstarSubset (readout)) { + psError(PSASTRO_ERR_DATA, false, "Can't determine an appropriate refstar subset\n"); + psFree (index); + psFree (view); + return false; + } + } + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; + + // activate all files except PSASTRO.OUTPUT + + pmFPAfileActivate (config->files, true, NULL); + pmFPAfileActivate (config->files, false, "PSASTRO.OUT.REFSTARS"); Index: /branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c (revision 17709) +++ /branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c (revision 17709) @@ -0,0 +1,99 @@ +# include "psphotInternal.h" + +bool psphotReadout(pmConfig *config, const pmFPAview *view) { + + psTimerStart ("psphotReadout"); + + // 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 (recipe, config, view)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } + + // 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"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + pmPSF *psf = psphotLoadPSF (config, view, recipe); + assert (psf); + + // load the source lists + // 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); + + // replace all sources + psphotRemoveSources (realSources, recipe); + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // add noise for real (subtracted) objects + psphotAddNoise (readout, realSources, recipe); + + // XXX fake sources should measure peak->x,y, force sources should not + psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); + ppSimDetections (significance, recipe, fakeSources); + ppSimDetections (significance, recipe, forceSources); + psFree (significance); + + // remove noise for subtracted objects (ie, return to normal noise level) + psphotSubNoise (readout, realSources, recipe); + + // replace all sources + psphotReplaceAll (realSources, recipe); + + // construct an initial model for each object + psphotGuessModels (readout, realSources, recipe, psf); + psphotGuessModels (readout, fakeSources, recipe, psf); + psphotGuessModels (readout, forceSources, recipe, psf); + + // linear fit to real + fake sources + psArray *sources = ppSimMergeSources (realSources, fakeSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAll (sources, recipe); + psFree (sources); // only frees the merged references + + // linear fit to real + forced sources + psArray *sources = ppSimMergeSources (realSources, forceSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAll (sources, recipe); + psFree (sources); // only frees the merged references + + // XXX do we need to measure aperture photometry corrections? + // XXX do we store the pre and post correction magnitudes? + // XXX for the fake sources, these must be identically zero + // XXX for the force sources, need the apresid to put the mags on the correct system + if (!psphotApResid (readout, forceSources, recipe, psf)) { + psLogMsg ("psphot", 3, "failed on psphotApResid"); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + } + + // calculate source magnitudes + pmReadout *background = psphotSelectBackground (config, view, false); + psphotMagnitudes(sources, recipe, psf, background); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (sources); + + // create the exported-metadata and free local data + return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); +}