Changeset 17935

Timestamp:

Jun 5, 2008, 3:29:09 AM (18 years ago)

Author:

eugene

Message:

added forced photometry code; compiles, not tested

Location:

branches/eam_branch_20080511/ppSim/src

Files:

• ppSimCreate.c (modified) (1 diff)
• ppSimInsertStars.c (modified) (2 diffs)
• ppSimLoadForceSources.c (modified) (1 diff)
• ppSimLoadSpots.c (modified) (5 diffs)
• ppSimLoop.c (modified) (1 diff)
• ppSimPhotomReadout.c (modified) (4 diffs)

branches/eam_branch_20080511/ppSim/src/ppSimCreate.c

@@ -29 +29 @@
   // 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
+  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.");

branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c

@@ -114 +114 @@
         source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
 
-        // XXX add the sources to a source array
-
         // insert the source flux in the image
         pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
         pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY);
-        psArrayAdd (sources, 100,source);
-        psFree(source);                 // Drop reference
 
         // Blow away the image parts of the source, which makes the memory explode
@@ -130 +126 @@
         RESET(source->psfFlux);
         RESET(source->blends);
+
+        // add the sources to the source array
+        psArrayAdd (sources, 100,source);
+        psFree(source);                 // Drop reference
     }
 

branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c

@@ -1 +1 @@
 # include "ppSim.h"
+
+// Reset a pointer: free and set to NULL
+#define RESET(PTR) \
+    psFree(PTR); \
+    PTR = NULL;
 
 psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view) {
 
-    psArray *result = psArrayAlloc(0);
-    return result;
+    bool status;
 
-# if (0)
+    pmFPAfile *outFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT");
+    psAssert(outFile, "missing PPSIM.OUTPUT");
+    psAssert(outFile->fpa, "missing fpr for PPSIM.OUTPUT");
 
-    // de-activate all files except PSASTRO.REFSTARS
-    pmFPAfileActivate (config->files, false, NULL);
-    pmFPAfileActivate (config->files, true, "PSASTRO.OUT.REFSTARS");
+    pmChip *outChip = pmFPAviewThisChip (view, outFile->fpa);
+    pmPSF *psf = psMetadataLookupPtr (&status, outChip->analysis, "PSPHOT.PSF");
+    assert (psf);
 
-    readout = pmFPAviewThisReadout (view, fpa);
-    if (! readout->data_exists) { continue; }
+    psArray *spots = psMetadataLookupPtr(&status, outFile->fpa->analysis, "FORCED.SPOTS");
+    PS_ASSERT_PTR_NON_NULL (spots, NULL);
 
-                // read WCS data from the corresponding header
-                pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
+    // in this program, we are looping over the output image, rather than the input as in ppImage
+    pmFPAfile *srcFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.REAL.SOURCES");
+    psAssert(srcFile, "missing PPSIM.REAL.SOURCES");
 
-                int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1");
-                int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2");
+    pmFPA     *fpa     = srcFile->fpa;
+    pmChip    *chip    = pmFPAviewThisChip (view, fpa);
+    // pmCell    *cell    = pmFPAviewThisCell (view, fpa);
+    pmReadout *readout = pmFPAviewThisReadout (view, fpa);
 
-                float minX = -fieldPadding*Nx;
-                float maxX = (1+fieldPadding)*Nx;
-                float minY = -fieldPadding*Ny;
-                float maxY = (1+fieldPadding)*Ny;
+    // XXX we should have only one cell and readout per chip, right?
 
-                // the refstars is a subset within range of this chip
-                psArray *refstars = psArrayAllocEmpty (100);
+    // read WCS from existing output from psphot & psastro
+    // XXX can we be more flexible?  ie, use the header of PPSIM.CHIP?
+    bool bilevelAstrometry = false;
+    status = psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
+    psAssert (status, "failed to set astrometry");
 
-                // 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]]);
+    // the pixel scale is only used to set the focal-plane coordinate system
+    // setting this value to 1 makes this system have the units of pixels
+    float pixelScale = 1.0;
+    status = psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry);
+    psAssert (status, "failed to set astrometry");
 
-                    psProject (ref->TP, ref->sky, fpa->toSky);
-                    psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
-                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
+    float minX = 0.0;
+    float maxX = readout->image->numCols;
+    float minY = 0.0;
+    float maxY = readout->image->numRows;
 
-                    // 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;
+    // the sources is a subset within range of this chip
+    psArray *sources = psArrayAllocEmpty (100);
 
-                    psArrayAdd (refstars, 100, ref);
-                skip:
-                    psFree (ref);
+    // Select the spots within range of this readout.  Project the spots to this chip
+    for (int i = 0; i < spots->n; i++) {
+        pmAstromObj *obj = spots->data[i];
 
-                    if (nMax && (refstars->n >= nMax)) break;
-                }
-                psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
+        psProject (obj->TP, obj->sky, fpa->toSky);
+        psPlaneTransformApply (obj->FP, fpa->fromTPA, obj->TP);
+        psPlaneTransformApply (obj->chip, chip->fromFPA, obj->FP);
 
-                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
-                psFree (refstars);
+        // limit the X,Y range of the objs to the selected chip
+        if (obj->chip->x < minX) continue;
+        if (obj->chip->x > maxX) continue;
+        if (obj->chip->y < minY) continue;
+        if (obj->chip->y > maxY) continue;
 
-                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;
+        // convert the pmAstromObj to pmSource
+
+        // instantiate a model for the PSF at this location, Io = 1.0
+        pmModel *model = pmModelFromPSFforXY (psf, obj->chip->x, obj->chip->y, 1.0);
+
+        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
+        // XXX use a fixed radius??
+        // float radius = model->modelRadius (model->params, roughNoise);
+        // radius = PS_MAX (radius, 1.0);
+        float radius = 5.0;
+
+        // construct a source, with model flux pixels set, based on the model
+        pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
+
+        // XXX set the mag & err values (should this be done in pmSourceFromModel?)
+        // XXX i should be applying the gain and the correct effective area
+        // psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
+        // psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+
+        // these are not really needed since we will be fitting for them
+        source->psfMag = NAN;
+        source->errMag = NAN;
+
+        // insert the source flux in the image
+        // XXX not sure the offset is really 0,0
+        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, 0.0, 0.0);
+        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, 0.0, 0.0);
+
+        // Blow away the image parts of the source, which makes the memory explode
+        RESET(source->pixels);
+        RESET(source->weight);
+        RESET(source->maskObj);
+        RESET(source->maskView);
+        RESET(source->modelFlux);
+        RESET(source->psfFlux);
+        RESET(source->blends);
+
+        psArrayAdd (sources, 100, source);
+        psFree(source);                 // Drop local reference
     }
-    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
+    psTrace ("psastro", 4, "Added %ld sources\n", sources->n);
 
-    // activate all files except PSASTRO.OUTPUT
-
-    pmFPAfileActivate (config->files, true, NULL);
-    pmFPAfileActivate (config->files, false, "PSASTRO.OUT.REFSTARS");
-# endif
-
-    return NULL;
+    return sources;
 }

branches/eam_branch_20080511/ppSim/src/ppSimLoadSpots.c

@@ -1 +1 @@
 # include "ppSim.h"
 
-// need to specify an alternative CATDIR for the force positions
-bool ppSimLoadSpots (psArray *spots, pmFPA *fpa, pmConfig *config) {
+// Load the relevant forced-photometry positions for this field.  This function does not determine
+// the pixel coordinates, merely the celestial coordinates for sources in the general vicinity.  The
+// sources are saved on the fpa->analysis metadata as FORCED.SPOTS.  We always create an entry; it
+// will be empty if no sources were selected or forced photometry is not desired.
+bool ppSimLoadSpots (pmFPA *fpa, pmConfig *config) {
 
     bool status;
+
+    psArray *spots = psArrayAllocEmpty (1);
+    psMetadataAddPtr (fpa->analysis, PS_LIST_TAIL, "FORCED.SPOTS", PS_META_REPLACE, "forced photometry positions", spots);
+    psFree (spots); // free the extra (local) reference; a copy remains on fpa->analysis
 
     psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe
@@ -15 +22 @@
     if (!astroRecipe) {
         psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE);
-        return NULL;
+        return false;
     }
 
@@ -22 +29 @@
     float scale   = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
 
+    char *catdir = psMetadataLookupStr(NULL, recipe, "FORCED.CATDIR");
+    char *photcode = psMetadataLookupStr(NULL, recipe, "FORCED.PHOTCODE");
+    float MAGmax = 100.0;
+
     // Size of FPA
     psRegion *bounds = ppSimFPABounds (fpa);
@@ -27 +38 @@
     psFree(bounds);
 
-    float x0fpa = 0.5*(bounds->x0 + bounds->x1);
-    float y0fpa = 0.5*(bounds->y0 + bounds->y1);
+    // modify the PSASTRO recipe to use the values desired for FORCED.PHOT.  we can use a view on
+    // the PSASTRO recipe because we are not calling psastro elsewhere in this program.  XXX need to
+    // use the WCS to define the overlap region
+    psMetadataAddStr(astroRecipe, PS_LIST_TAIL, "DVO.CATDIR",  PS_META_REPLACE, "", catdir);
+    psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", ra0 - radius);
+    psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "RA_MAX",  PS_META_REPLACE, "", ra0 + radius);
+    psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", dec0 - radius);
+    psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", dec0 + radius);
+    
+    // XXX these two values are somewhat bogus: can we have psastroLoadRefstars optionally ignore photcode and MAGmax?
+    psMetadataAddStr(astroRecipe, PS_LIST_TAIL, "DVO.GETSTAR.PHOTCODE",  PS_META_REPLACE, "", photcode);
+    psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DVO.GETSTAR.MAG.MAX",  PS_META_REPLACE, "", MAGmax);
 
-    // XXX need to use the WCS to define the overlap region
-    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius);
-    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 + radius);
-    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", dec0 - radius);
-    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius);
     psArray *refStars = psastroLoadRefstars(config);
     if (!refStars || refStars->n == 0) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars.");
         psFree(refStars);
-        return NULL;
+        return false;
     }
     psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n);
 
     // convert the pmAstromObj sources to the storage format used by the CMF files
-
-    long oldSize = spots->n;
-    spots = psArrayRealloc (spots, oldSize + refStars->n);
-
-    // Conversion loop
+    psArrayRealloc (spots, refStars->n);
     for (long i = 0; i < refStars->n; i++) {
         ppSimStar *spot = ppSimStarAlloc ();
@@ -62 +74 @@
         // spot->y = sin(pa) * xi + cos(pa) * eta + y0fpa;
 
-        spots->data[oldSize + i] = star;
+        spots->data[i] = spot;
 
         psTrace("ppSim", 10, "Adding catalogue star: %.1f,%.1f --> %.2f\n", spot->x, spot->y, ref->Mag);
     }
 
-    // XXX these need to be saved on PSPHOT.INPUT.CMF
     return true;
 }

branches/eam_branch_20080511/ppSim/src/ppSimLoop.c

@@ -27 +27 @@
     int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y
 
-    psArray *spots = psArrayAllocEmpty (1);
     psArray *stars = psArrayAllocEmpty (1);
     psArray *galaxies = psArrayAllocEmpty (1);
     if (type == PPSIM_TYPE_OBJECT) {
-        // Load forced-photometry positions
-        if (!ppSimLoadSpots (spots, fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load forced-photometry spots");
+        // Load forced-photometry positions (these are placed on fpa->analysis for use in ppSimPhotomReadout)
+        if (!ppSimLoadSpots (fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load forced-photometry spots");
 
         // Load catalogue stars

branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c

@@ -37 +37 @@
     }
 
+# if 0    
     // set the photcode for this image
     if (!psphotAddPhotcode (recipe, config, view, "PPSIM.CHIP")) {
@@ -42 +43 @@
         return false;
     }
+# endif
 
     // find the currently selected readout. 
@@ -119 +121 @@
     psFree (sources); // only frees the merged references
 
-    // linear fit to real + forced sources
-    sources = ppSimMergeSources (realSources, forceSources);
-    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract
-    psphotReplaceAllSources (sources, recipe);
-    psFree (sources); // only frees the merged references
+    // linear fit to real + forced sources (if any were selected)
+    if (forceSources->n) {
+      sources = ppSimMergeSources (realSources, forceSources);
+      psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract
+      psphotReplaceAllSources (sources, recipe);
+      psFree (sources); // only frees the merged references
+    }
 
     // XXX do we need to measure aperture photometry corrections?
@@ -157 +161 @@
     psMetadataAdd (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "fake photometry ", fakeSources);
 
-    // psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PPSIM.FORCE.SOURCES", PS_DATA_ARRAY, "forced photometry ", forceSources);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PPSIM.FORCE.SOURCES", PS_DATA_ARRAY, "forced photometry ", forceSources);
     psFree (forceSources);
     return true;