Changeset 14667

Timestamp:

Aug 25, 2007, 11:57:36 AM (19 years ago)

Author:

eugene

Message:

added galaxies, some tweaks

Location:

trunk/ppSim/src

Files:

• Makefile.am (modified) (2 diffs)
• ppSim.h (modified) (3 diffs)
• ppSimArguments.c (modified) (2 diffs)
• ppSimCreate.c (modified) (2 diffs)
• ppSimInsertGalaxies.c (added)
• ppSimInsertSources.c (modified) (2 diffs)
• ppSimLoadStars.c (modified) (3 diffs)
• ppSimLoop.c (modified) (7 diffs)
• ppSimMakeBias.c (modified) (1 diff)
• ppSimMakeGalaxies.c (added)
• ppSimMakeStars.c (modified) (2 diffs)
• ppSimStars.c (modified) (2 diffs)

trunk/ppSim/src/Makefile.am

@@ -9 +9 @@
         ppSimLoadStars.c        \
         ppSimMakeStars.c        \
+        ppSimMakeGalaxies.c     \
         ppSimMakeBiassec.c      \
         ppSimMakeBias.c         \
@@ -14 +15 @@
         ppSimMakeSky.c          \
         ppSimInsertSources.c    \
+        ppSimInsertGalaxies.c   \
         ppSimAddOverscan.c      \
         ppSimAddNoise.c         \

trunk/ppSim/src/ppSim.h

@@ -48 +48 @@
 } ppSimStar;
 
+typedef struct {
+    double ra;
+    double dec;
+    float mag;
+    float x;
+    float y;
+    float flux;
+
+    float peak;
+    float Rmaj;
+    float Rmin;
+    float theta;
+    float index;
+} ppSimGalaxy;
 
 ppSimStar *ppSimStarAlloc ();
+ppSimGalaxy *ppSimGalaxyAlloc ();
 
 /// Parse command-line arguments
@@ -85 +100 @@
                      );
 
-psArray *ppSimLoadStars (pmFPA *fpa, pmConfig *config);
+bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config);
 bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng);
 
@@ -117 +132 @@
 bool ppSimSetPSF (pmChip *chip, pmConfig *config);
 
+bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng);
+bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config);
+
 #endif

trunk/ppSim/src/ppSimArguments.c

@@ -71 +71 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist");
 
-    if (!config->camera) {
-        psErrorStackPrint(stderr, "A camera name must be specified using the -camera option.");
+    // only one of -camera and -file is needed
+    bool status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
+    if (!config->camera && !status) {
+        psErrorStackPrint(stderr, "A camera name (-camera NAME) or an image (-file NAME) must be specified");
         usage(argv[0], arguments, config);
     }
@@ -82 +84 @@
     psString formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format
     if (formatName) {
+        // XXX delay the config below until ppSimCreate?
         config->formatName = psMemIncrRefCounter(formatName);
 

trunk/ppSim/src/ppSimCreate.c

@@ -9 +9 @@
 pmFPAfile *ppSimCreate(pmConfig *config)
 {
+    bool status;
+    bool simImage = false;
     PS_ASSERT_PTR_NON_NULL(config, NULL);
+    pmFPA *fpa = NULL;
 
-    pmFPA *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;
+    // 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) {
+        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 {
+        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;
     }
 
@@ -56 +72 @@
     }
     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

trunk/ppSim/src/ppSimInsertSources.c

@@ -5 +5 @@
     bool mdok;
 
+    assert (readout);
     assert (stars);
-    assert (readout);
-    
+
+    if (!stars->n) { return true; }
+
     // XXX is this needed?
     // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources
@@ -88 +90 @@
         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;
+
+        source->psfMag = -2.5*log10(star->flux);
+        source->errMag = sqrt(Area*PS_SQR(roughNoise) + star->flux) / star->flux;       
+        
         // XXX add the sources to a source array
 

trunk/ppSim/src/ppSimLoadStars.c

@@ -1 +1 @@
 #include "ppSim.h"
 
-psArray *ppSimLoadStars (pmFPA *fpa, pmConfig *config) {
+bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config) {
+
+    bool mdok;
+    assert (stars);
+
+    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
+
+    bool starsReal = psMetadataLookupBool(&mdok, recipe, "STARS.REAL"); // Density of fakes
+    if (!starsReal) return true;
 
     // Read catalogue stars using psastro
@@ -38 +46 @@
     psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n);
 
-    psArray *stars = psArrayAlloc (refStars->n);
+    long oldSize = stars->n;
+    stars = psArrayRealloc (stars, refStars->n);
 
     // Conversion loop
@@ -61 +70 @@
         star->flux = powf(10.0, -0.4 * (star->mag - zp)) * expTime;
         star->peak = star->flux / sqrt(2.0*M_PI) / seeing;
-        stars->data[i] = star;
+        stars->data[oldSize + i] = star;
     }
     psFree(refStars);

trunk/ppSim/src/ppSimLoop.c

@@ -17 +17 @@
 
     // Load catalogue stars
-    psArray *stars = NULL;
+    psArray *stars = psArrayAllocEmpty (1);
     if (type == PPSIM_TYPE_OBJECT) {
-        stars = ppSimLoadStars (fpa, config);
+        ppSimLoadStars (stars, fpa, config);
     }
 
@@ -25 +25 @@
     if (type == PPSIM_TYPE_OBJECT) {
         ppSimMakeStars (stars, fpa, config, rng);
+    }
+
+    // Add random galaxies
+    psArray *galaxies = psArrayAllocEmpty (1);
+    if (type == PPSIM_TYPE_OBJECT) {
+        ppSimMakeGalaxies (galaxies, fpa, config, rng);
     }
 
@@ -62 +68 @@
         while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
 
-            // Size, position and orientation of cell
-            int numCols = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE") / binning;
-            int numRows = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE") / binning;
-
+            // check that we are able to work with this cell (readdir must be 1)
             int readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); // Read direction
             if (readdir != 1) {
@@ -73 +76 @@
             }
 
-            psVector *biasCols = ppSimMakeBiassec (cell, config);
+            // if no readout exists, generate a single readout (upgrade to allow multiple
+            // readouts?)
+            if (cell->readouts->n == 0) {
+                // Size, position and orientation of cell
+                int numCols = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE") / binning;
+                int numRows = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE") / binning;
 
-            int numReadouts = 1;
-            for (int i = 0; i < numReadouts; i++) {
+                // generate a new readout for this cell
                 pmReadout *readout = pmReadoutAlloc(cell); // Readout within cell
 
@@ -83 +90 @@
                 readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
 
+                psImageInit (readout->image, 0.0);
+                psImageInit (readout->weight, 0.0);
+            }
+
+            psVector *biasCols = ppSimMakeBiassec (cell, config);
+
+            for (int i = 0; i < cell->readouts->n; i++) {
+
+                pmReadout *readout = cell->readouts->data[i];
+                assert (readout);
+
+                // if we have not read in a weight or generated a fake image above, we need to
+                // build one here
+                if (!readout->weight) {
+                    if (!pmReadoutGenerateWeight(readout, true)) {
+                        psError (PS_ERR_UNKNOWN, false, "trouble creating weight");
+                        return false;
+                    }
+                }
+
                 psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects
                 if (type == PPSIM_TYPE_OBJECT) {
-                    expCorr = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+                    expCorr = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
                 }
 
                 psVector *biasRows = ppSimMakeBias (readout, config, rng);
                 if (type == PPSIM_TYPE_BIAS) goto done;
-                
+        
                 ppSimMakeDark (readout, config);
                 if (type == PPSIM_TYPE_DARK) goto done;
@@ -99 +126 @@
                 if (type == PPSIM_TYPE_OBJECT) {
                     ppSimInsertSources (readout, expCorr, stars, config);
+                }
+
+                if (type == PPSIM_TYPE_OBJECT) {
+                    ppSimInsertGalaxies (readout, expCorr, galaxies, config);
                 }
 
@@ -111 +142 @@
                 readout->parent->data_exists = true;
                 readout->parent->parent->data_exists = true;
-                psFree(readout);        // Drop reference
             }
             psFree(biasCols);

trunk/ppSim/src/ppSimMakeBias.c

@@ -46 +46 @@
 
             // Bias level
-            signal->data.F32[y][x] = biasRows->data.F32[y];
-            variance->data.F32[y][x] = PS_SQR(readnoise);
+            signal->data.F32[y][x] += biasRows->data.F32[y];
+            variance->data.F32[y][x] += PS_SQR(readnoise);
 
         }

trunk/ppSim/src/ppSimMakeStars.c

@@ -3 +3 @@
 bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
 
+    bool mdok;
     assert (stars);
 
-    bool mdok;
+    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
 
-    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
+    bool starsFake = psMetadataLookupBool(&mdok, recipe, "STARS.FAKE"); // Density of fakes
+    if (!starsFake) return true;
 
     float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
@@ -45 +47 @@
         psLogMsg("ppSim", PS_LOG_INFO,
                  "Image noise is above brightest random star --- no random stars added.");
-    } else {
+        return true;
+    }
 
-        // Normalisation, set by the specified stellar density at the specified bright magnitude
-        float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
-            powf(bright, starsLum);
+    // Normalisation, set by the specified stellar density at the specified bright magnitude
+    float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
+        powf(bright, starsLum);
 
-        // Total number of stars down to the faint flux end
-        long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
+    // Total number of stars down to the faint flux end
+    long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
 
-        psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
-                 -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
+    psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
+             -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
 
-        long oldSize = stars->n;
-        psArrayRealloc (stars, stars->n + num);
+    long oldSize = stars->n;
+    psArrayRealloc (stars, stars->n + num);
 
-        for (long i = 0; i < num; i++) {
-            ppSimStar *star = ppSimStarAlloc ();
+    for (long i = 0; i < num; i++) {
+        ppSimStar *star = ppSimStarAlloc ();
 
-            // make fpa center of distribution
-            star->x    = psRandomUniform(rng) * xSize; // x position
-            star->y    = psRandomUniform(rng) * ySize; // y position
-            star->peak = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
-            star->flux = star->peak * sqrt(2.0*M_PI)*seeing;
-            stars->data[oldSize + i] = star;
-        }
-        stars->n = stars->nalloc;
+        // make fpa center of distribution
+        star->x    = psRandomUniform(rng) * xSize; // x position
+        star->y    = psRandomUniform(rng) * ySize; // y position
+
+        star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
+        star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));
+
+        stars->data[oldSize + i] = star;
     }
+    stars->n = stars->nalloc;
+
     return true;
 }

trunk/ppSim/src/ppSimStars.c

@@ -1 +1 @@
 # include "ppSim.h"
-
-#define STAR_RANGE_MIN 4.5              // Minimum pixel range for adding stars, sigma
-#define STAR_BG_FRAC 0.1                // Fraction of background error to go out to when adding stars
-
-// Return star flux at a position
-static inline float starFlux(float x, float y, // Position of interest
-                             float x0, float y0, // Centre of star
-                             float seeing // Seeing FWHM (pixels)
-    )
-{
-#ifdef STAR_GAUSSIAN
-    // Gaussian star
-    return exp(-0.5 * (PS_SQR(x - x0) + PS_SQR(y - y0)) / PS_SQR(seeing));
-#else
-    // Waussian star
-    float z = (PS_SQR(x - x0) + PS_SQR(y - y0)) / (2.0 * PS_SQR(seeing));
-    return 1.0 / (1.0 + z + PS_SQR(z));
-#endif
-}
-
-// Return size of star in pixels
-static inline int starSize(float noise, float peak, float seeing)
-{
-#ifdef STAR_GAUSSIAN
-    // Gaussian star (solving Gaussian)
-    float target = STAR_BG_FRAC * seeing * sqrt(M_2_PI) * noise / peak;
-    return target < 1.0 ? 2.0 * seeing * sqrtf(-logf(target)) + 0.5 : seeing * STAR_RANGE_MIN;
-#else
-    // Waussian star (solving Waussian where z >> 1 and peak/sqrt(bg) >> 1)
-    float target = STAR_BG_FRAC * noise / peak;
-    return PS_MAX(sqrtf(1.0 / target) + 0.5, seeing * STAR_RANGE_MIN);
-#endif
-}
-
-// Add a star into the signal and variance images
-bool ppSimInsertStar(psImage *signal,       // Signal image, to which to add star
-                     psImage *variance,     // Variance image, to which to add star
-                     float x0, float y0,    // Position of star
-                     float peak,            // Peak flux of star
-                     float noise,           // Rough noise estimate
-                     float seeing,          // Seeing for star
-                     const psImage *correction // Exposure correction as a function of position
-    )
-{
-    int size = starSize(noise, peak, seeing); // Size of star
-
-    // Range in image pixels on which to add star
-    int xMin = PS_MAX(0, x0 - size);
-    int xMax = PS_MIN(signal->numCols - 1, x0 + size);
-    int yMin = PS_MAX(0, y0 - size);
-    int yMax = PS_MIN(signal->numRows - 1, y0 + size);
-
-    for (int y = yMin; y <= yMax; y++) {
-        for (int x = xMin; x <= xMax; x++) {
-            float star = peak * correction->data.F32[y][x] * starFlux(x, y, x0, y0, seeing); // Star flux
-            signal->data.F32[y][x] += star;
-            variance->data.F32[y][x] += star;
-        }
-    }
-    return true;
-}
 
 void ppSimStarFree(ppSimStar *star)
@@ -74 +13 @@
     return star;
 }
+
+void ppSimGalaxyFree(ppSimGalaxy *galaxy)
+{
+    return;
+}
+
+ppSimGalaxy *ppSimGalaxyAlloc () {
+
+    ppSimGalaxy *galaxy = (ppSimGalaxy *) psAlloc(sizeof(ppSimGalaxy));
+    psMemSetDeallocator(galaxy, (psFreeFunc) ppSimGalaxyFree);
+
+    return galaxy;
+}