Changeset 25383

Timestamp:

Sep 15, 2009, 12:45:01 PM (17 years ago)

Author:

Paul Price

Message:

Merging branches/pap (detection efficiency development) into trunk.

Location:

trunk

Files:

• . (modified) (1 prop)
• archive/ducttape (modified) (1 prop)
• catalyst (modified) (1 prop)
• catalyst/2007.0-specs (modified) (1 prop)
• catalyst/2008.0-specs (modified) (1 prop)
• catalyst/2008.0-specs/hardened (modified) (1 prop)
• doc/misc/docgen.pl (modified) (1 prop)
• hardware (modified) (1 prop)
• ippScripts/scripts/magic_destreak_revert.pl (modified) (1 prop)
• ippTests/tap (modified) (1 prop)
• ippconfig/gpc1/camera.config (modified) (1 diff)
• ippconfig/recipes/psphot.config (modified) (1 diff)
• psLib/src/fits/psFits.c (modified) (4 diffs)
• psLib/src/fits/psFits.h (modified) (1 diff)
• psLib/src/imageops/psImageConvolve.c (modified) (4 diffs)
• psLib/src/imageops/psImageConvolve.h (modified) (1 diff)
• psLib/src/types/psMetadata.c (modified) (4 diffs)
• psLib/src/types/psMetadata.h (modified) (1 diff)
• psModules/src/camera/pmFPAMosaic.c (modified) (13 diffs)
• psModules/src/camera/pmReadoutFake.c (modified) (2 diffs, 1 prop)
• psModules/src/camera/pmReadoutFake.h (modified) (1 diff)
• psModules/src/objects/Makefile.am (modified) (2 diffs)
• psModules/src/objects/pmDetEff.c (copied) (copied from branches/pap/psModules/src/objects/pmDetEff.c )
• psModules/src/objects/pmDetEff.h (copied) (copied from branches/pap/psModules/src/objects/pmDetEff.h )
• psModules/src/objects/pmSourceIO.c (modified) (8 diffs)
• psModules/src/psmodules.h (modified) (1 diff)
• psphot/doc/efficiency.txt (modified) (1 prop)
• psphot/src/Makefile.am (modified) (1 diff)
• psphot/src/psphot.h (modified) (4 diffs)
• psphot/src/psphotAddNoise.c (modified) (3 diffs)
• psphot/src/psphotEfficiency.c (copied) (copied from branches/pap/psphot/src/psphotEfficiency.c )
• psphot/src/psphotFake.c (modified) (1 prop)
• psphot/src/psphotFitSourcesLinear.c (modified) (1 diff)
• psphot/src/psphotMagnitudes.c (modified) (2 diffs)
• psphot/src/psphotReadout.c (modified) (2 diffs)
• psphot/src/psphotReadoutMinimal.c (modified) (1 diff)
• psphot/src/psphotReplaceUnfit.c (modified) (1 diff)
• psphot/src/psphotSignificanceImage.c (modified) (1 diff)

trunk/ippconfig/gpc1/camera.config

@@ -132 +132 @@
   CMF.XSRC STR {CHIP.NAME}.xsrc # use .PSF and .EXT?
   CMF.XFIT STR {CHIP.NAME}.xfit # use .PSF and .EXT?
+  CMF.DETEFF STR {CHIP.NAME}.deteff
 
   PSF.HEAD  STR {CHIP.NAME}.hdr

trunk/ippconfig/recipes/psphot.config

@@ -265 +265 @@
 PSPHOT.CR.NSIGMA.SOFTEN             F32   0.0025          # Softening parameter for weights
 
+# Detection efficiency
+EFF.NUM                             S32   500                   # Number of fake sources per bin
+@EFF.MAG                            F32   -2.0 -1.0 -0.5 -0.25 -0.1 -0.05 0.0 0.05 0.1 0.25 0.5 1.0 2.0 # Magnitude of fake sources relative to limit
 
 # Recipe overrides for CHIP

trunk/psLib/src/fits/psFits.c

@@ -344 +344 @@
 // Therefore, we implement our own version of moving to an extension specified by name.  The pure cfitsio
 // version is used if "conventions.compression" handling is turned off in the psFits structure.
-bool psFitsMoveExtName(const psFits* fits,
-                       const char* extname)
+static bool fitsMoveExtName(const psFits* fits, // FITS file
+                            const char* extname, // Extension name
+                            bool errors // Generate errors?
+    )
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
@@ -356 +358 @@
         // User wants to use cfitsio.  Good luck to them!
         if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
-            psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            }
             return false;
         }
@@ -378 +382 @@
         if (fits_movabs_hdu(fits->fd, i, &hdutype, &status)) {
             // We've run off the end
-            psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            }
             return false;
         }
@@ -406 +412 @@
     }
     psAbort("Should never reach here.");
+}
+
+
+bool psFitsMoveExtName(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, true);
+}
+
+bool psFitsMoveExtNameClean(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, false);
 }
 

trunk/psLib/src/fits/psFits.h

@@ -245 +245 @@
 );
 
+/** Moves the FITS HDU to the specified extension name without generating errors.
+ *
+ *  @return bool        TRUE if the extension name was found and move was
+ *                      successful, otherwise FALSE
+ */
+bool psFitsMoveExtNameClean(
+    const psFits* fits,                ///< the psFits object to move
+    const char* extname                ///< the extension name
+);
+
 /** Moves the FITS HDU to the specified extension number
  *

trunk/psLib/src/imageops/psImageConvolve.c

@@ -678 +678 @@
 }
 
+static bool imageSmoothMaskPixels(psVector *out, const psImage *image, const psImage *mask,
+                                  psImageMaskType maskVal, const psVector *x, const psVector *y,
+                                  const psVector *gaussNorm, float minGauss, int size, int start, int stop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int xLast = numCols - 1, yLast = numRows - 1; // Last index
+    for (int i = start; i < stop; i++) {
+        int xPix = x->data.S32[i], yPix = y->data.S32[i]; // Pixel coordinates for smoothing
+
+        int yMin = PS_MAX(yPix - size, 0);
+        int yMax = PS_MIN(yPix + size, yLast);
+        int xMin = PS_MAX(xPix - size, 0);
+        int xMax = PS_MIN(xPix + size, xLast);
+
+        const float *yGauss = &gauss[yMin - yPix];
+
+        double ySumIG = 0.0, ySumG = 0.0;
+        for (int v = yMin; v <= yMax; v++, yGauss++) {
+            const float *xGauss = &gauss[xMin - xPix];
+            double xSumIG = 0.0, xSumG = 0.0;
+            const psImageMaskType *maskData = &mask->data.PS_TYPE_IMAGE_MASK_DATA[v][xMin];
+            const psF32 *imageData = &image->data.F32[v][xMin];
+            for (int u = xMin; u <= xMax; u++, xGauss++, imageData++, maskData++) {
+                if (*maskData & maskVal) {
+                    continue;
+                }
+                xSumIG += *imageData * *xGauss;
+                xSumG += *xGauss;
+            }
+            if (xSumG > minGauss) {
+                ySumIG += xSumIG * *yGauss;
+                ySumG += xSumG * *yGauss;
+            }
+        }
+
+        out->data.F32[i] = ySumG > minGauss ? ySumIG / ySumG : NAN;
+    }
+
+    return true;
+}
+
+static bool psImageSmoothMaskPixelsThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "programming error: no job arguments");
+    psAssert(job->args->n == 11, "programming error: wrong number of job arguments");
+
+    psVector *out = job->args->data[0]; // Output vector
+    const psImage *image  = job->args->data[1]; // Input image
+    const psImage *mask   = job->args->data[2]; // Input mask
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3], PS_TYPE_IMAGE_MASK_DATA);
+    const psVector *x = job->args->data[4];
+    const psVector *y = job->args->data[5];
+    const psVector *gaussNorm = job->args->data[6];
+    float minGauss = PS_SCALAR_VALUE(job->args->data[7], F32);
+    int size = PS_SCALAR_VALUE(job->args->data[8], S32);
+    int start = PS_SCALAR_VALUE(job->args->data[9], S32);
+    int stop = PS_SCALAR_VALUE(job->args->data[10], S32);
+    return imageSmoothMaskPixels(out, image, mask, maskVal, x, y, gaussNorm,
+                                 minGauss, size, start, stop);
+}
+
+
+psVector *psImageSmoothMaskPixels(const psImage *image, const psImage *mask, psImageMaskType maskVal,
+                                  const psVector *x, const psVector *y,
+                                  float sigma, float numSigma, float minGauss)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_IMAGE_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+
+    int num = x->n;                     // Number of pixels to smooth
+    psVector *out = psVectorAlloc(num, PS_TYPE_F32); // Output results
+
+    // Generate normalized gaussian
+    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
+
+    // Columns
+    if (threaded) {
+        int numThreads = psThreadPoolSize(); // Number of threads
+        int delta = (numThreads) ? num / numThreads + 1 : num; // Block of cols to do at once
+        for (int start = 0; start < num; start += delta) {
+            int stop = PS_MIN(start + delta, num);  // End of block
+
+            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
+            psArrayAdd(job->args, 1, out);
+            psArrayAdd(job->args, 1, (psImage*)image);
+            psArrayAdd(job->args, 1, (psImage*)mask);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            psArrayAdd(job->args, 1, (psVector*)x);
+            psArrayAdd(job->args, 1, (psVector*)y);
+            psArrayAdd(job->args, 1, gaussNorm);
+            PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(gaussNorm);
+                psFree(out);
+                return NULL;
+            }
+            psFree(job);
+        }
+    } else if (!imageSmoothMaskPixels(out, image, mask, maskVal, x, y,
+                                      gaussNorm, minGauss, size, 0, num)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to smooth pixels.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    psFree(gaussNorm);
+    return out;
+}
+
+
 // Smooth an image with masked pixels
 // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
@@ -1325 +1458 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaRows = (numThreads) ? numRows / numThreads : numRows; // Block of rows to do at once
+        int deltaRows = (numThreads) ? numRows / numThreads + 1 : numRows; // Block of rows to do at once
         for (int start = 0; start < numRows; start += deltaRows) {
-            int stop = PS_MIN (start + deltaRows, numRows);  // end of row block
+            int stop = PS_MIN(start + deltaRows, numRows);  // end of row block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1363 +1496 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaCols = (numThreads) ? numCols / numThreads : numCols; // Block of cols to do at once
+        int deltaCols = (numThreads) ? numCols / numThreads + 1 : numCols; // Block of cols to do at once
         for (int start = 0; start < numCols; start += deltaCols) {
-            int stop = PS_MIN (start + deltaCols, numCols);  // end of col block
+            int stop = PS_MIN(start + deltaCols, numCols);  // end of col block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1486 +1619 @@
             psFree(task);
         }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9);
+            task->function = &psImageSmoothMaskPixelsThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
     } else if (!set && threaded) {
         psThreadTaskRemove("PSLIB_IMAGE_CONVOLVE_MASK");

trunk/psLib/src/imageops/psImageConvolve.h

@@ -222 +222 @@
     );
 
+/// Smooth particular pixels on an image, allowing for masked pixels
+///
+/// Applies a circularly symmetric Gaussian smoothing first in x and then in y
+/// directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
+psVector *psImageSmoothMaskPixels(
+    const psImage *image,               ///< Input image (F32)
+    const psImage *mask,                ///< Mask image
+    psImageMaskType maskVal,            ///< Value to mask
+    const psVector *x,                  ///< x coordinates
+    const psVector *y,                  ///< y coordinates
+    float sigma,                        ///< Width of the smoothing kernel (pixels)
+    float numSigma,                     ///< Size of the smoothing box (sigma)
+    float minGauss                      ///< Minimum fraction of Gaussian to accept
+    );
 
 psImage *psImageSmoothMask_Threaded(psImage *output,

trunk/psLib/src/types/psMetadata.c

@@ -996 +996 @@
     if (metadataItem->type == PS_DATA_METADATA_MULTI) {
         // if multiple keys found, use the first.
-        //        metadataItem = (psMetadataItem*)((metadataItem->data.list)->head);
         metadataItem = (psMetadataItem*)(metadataItem->data.list->head->data);
         if (status) {
@@ -1082 +1081 @@
     } \
     \
-    /* psFree(metadataItem); currently, the lookup doesn't increment the ref count */ \
     return value; \
 }
@@ -1101 +1099 @@
 psMetadataLookupNumTYPE(VectorMaskType,VectorMask)
 psMetadataLookupNumTYPE(ImageMaskType,ImageMask)
+
+#define psMetadataLookupPtrTYPE(TYPENAME,NAME,TYPE,VAL) \
+TYPENAME psMetadataLookup##NAME(bool *status, const psMetadata *md, const char *key) \
+{ \
+    PS_ASSERT_METADATA_NON_NULL(md, NULL); \
+    PS_ASSERT_STRING_NON_EMPTY(key, NULL); \
+    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); /* The item of interest */ \
+    if (!item) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (item->type == PS_DATA_METADATA_MULTI) { \
+        /* if multiple keys found, use the first. */ \
+        item = item->data.list->head->data; \
+    } \
+    if (item->type != TYPE) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (status) { \
+        *status = true; \
+    } \
+    return item->data.VAL; \
+}
+
+psMetadataLookupPtrTYPE(psMetadata*, Metadata, PS_DATA_METADATA, md)
+psMetadataLookupPtrTYPE(psString, Str, PS_DATA_STRING, str)
+psMetadataLookupPtrTYPE(psTime*, Time, PS_DATA_TIME, V)
+psMetadataLookupPtrTYPE(psVector*, Vector, PS_DATA_VECTOR, V)
+
 
 psMetadataItem* psMetadataGet(const psMetadata *md,
@@ -1257 +1295 @@
 }
 
-psMetadata *psMetadataLookupMetadata(bool *status,
-                                     const psMetadata *md,
-                                     const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    psMetadata *value = NULL;  // The value to return
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_METADATA) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.md; // The requested metadata
-    }
-    return value;
-}
-
-psTime *psMetadataLookupTime(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key);
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-        return NULL;
-    }
-
-    if (item->type != PS_DATA_TIME) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_TIME, as expected.\n", key);
-        }
-        return NULL;
-    }
-
-    // We have the requested metadata
-    if (status) {
-        *status = true;
-    }
-
-    return item->data.V;
-}
-
-
-psString psMetadataLookupStr(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    //    char *value = NULL;   // The value to return
-    psString value = NULL;
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_STRING) {
-        // The value at the key isn't of the desired type
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_STRING, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.V;
-    }
-    return value;
-}
+
+
+
 
 psList *psMetadataKeys(psMetadata *md)

trunk/psLib/src/types/psMetadata.h

@@ -717 +717 @@
 PS_METADATA_LOOKUP_TYPE_DECL(Ptr, psPtr);
 PS_METADATA_LOOKUP_TYPE_DECL(Str, psString);
+PS_METADATA_LOOKUP_TYPE_DECL(Vector, psVector*);
 PS_METADATA_LOOKUP_TYPE_DECL(Metadata, psMetadata*);
 PS_METADATA_LOOKUP_TYPE_DECL(Time, psTime*);

trunk/psModules/src/camera/pmFPAMosaic.c

@@ -740 +740 @@
 static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                        psImage **mosaicMask, // The mosaic mask, to be returned
-                       psImage **mosaicVariances, // The mosaic variances, to be returned
+                       psImage **mosaicVariance, // The mosaic variance, to be returned
                        int *xBinChip, int *yBinChip, // The binning in x and y, to be returned
                        const pmChip *chip, // Chip to mosaic
@@ -749 +749 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicVariances);
+    assert(mosaicVariance);
     assert(xBinChip);
     assert(yBinChip);
@@ -826 +826 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
-        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
+        *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank);
     }
@@ -847 +847 @@
 static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                       psImage **mosaicMask, // The mosaic mask, to be returned
-                      psImage **mosaicVariances, // The mosaic variances, to be returned
+                      psImage **mosaicVariance, // The mosaic variance, to be returned
                       int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned
                       const pmFPA *fpa,  // FPA to mosaic
@@ -857 +857 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicVariances);
+    assert(mosaicVariance);
     assert(xBinFPA);
     assert(yBinFPA);
@@ -960 +960 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
-        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
+        *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank);
     }
@@ -1025 +1025 @@
     psImage *mosaicImage   = NULL;      // The mosaic image
     psImage *mosaicMask    = NULL;      // The mosaic mask
-    psImage *mosaicVariances = NULL;      // The mosaic variances
+    psImage *mosaicVariance = NULL;      // The mosaic variances
 
     // Find the HDU
@@ -1052 +1052 @@
         }
         if (hdu->variances) {
-            mosaicVariances = psImageSubset(hdu->variances->data[0], bounds);
-            if (!mosaicVariances) {
+            mosaicVariance = psImageSubset(hdu->variances->data[0], bounds);
+            if (!mosaicVariance) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to select variance pixels.\n");
                 return false;
@@ -1061 +1061 @@
         // Case 2 --- we need to mosaic by cut and paste
         psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n");
-        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) {
+        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariance, &xBin, &yBin, source, targetCell, blank)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n");
             return false;
@@ -1069 +1069 @@
     }
     psTrace("psModules.camera", 1, "xBin,yBin: %d,%d\n", xBin, yBin);
+
 
     // Set the concepts for the target cell
@@ -1090 +1091 @@
     target->parent->concepts = psMetadataCopy(target->parent->concepts, source->parent->concepts); // FPA lvl
 
+    // Average the covariances
+    psList *covariances = psListAlloc(NULL); // Input covariance matrices
+    for (int i = 0; i < source->cells->n; i++) {
+        pmCell *cell = source->cells->data[i]; // Cell of interest
+        if (!cell || !cell->data_exists) {
+            continue;
+        }
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+        if (!ro || !ro->covariance) {
+            continue;
+        }
+        psListAdd(covariances, PS_LIST_TAIL, ro->covariance);
+    }
+    psKernel *mosaicCovariance = NULL;  // Covariance for mosaic
+    if (psListLength(covariances) > 0) {
+        psArray *covarArray = psListToArray(covariances); // Array with covariances
+        mosaicCovariance = psImageCovarianceAverage(covarArray);
+        psFree(covarArray);
+    }
+    psFree(covariances);
+
     // Now make a new readout to go in the target cell
     pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
     newReadout->image  = mosaicImage;
     newReadout->mask   = mosaicMask;
-    newReadout->variance = mosaicVariances;
+    newReadout->variance = mosaicVariance;
+    newReadout->covariance = mosaicCovariance;
     psFree(newReadout);                 // Drop reference
 
@@ -1334 +1357 @@
     target->concepts = psMetadataCopy(target->concepts, source->concepts);
 
+    // Average the covariances
+    psList *covariances = psListAlloc(NULL); // Input covariance matrices
+    for (int i = 0; i < covariances->n; i++) {
+        pmChip *chip = chips->data[i];  // Chip of interest
+        if (!chip || !chip->data_exists) {
+            continue;
+        }
+        psArray *cells = chip->cells;   // Cells in chip
+        for (long j = 0; j < cells->n; j++) {
+            pmCell *cell = cells->data[i]; // Cell of interest
+            if (!cell || !cell->data_exists) {
+                continue;
+            }
+            pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+            if (!ro || !ro->covariance) {
+                continue;
+            }
+            psListAdd(covariances, PS_LIST_TAIL, ro->covariance);
+        }
+    }
+    psKernel *mosaicCovariances = NULL; // Covariance for mosaic
+    if (psListLength(covariances) > 0) {
+        psArray *covarArray = psListToArray(covariances); // Array with covariances
+        mosaicCovariances = psImageCovarianceAverage(covarArray);
+        psFree(covarArray);
+    }
+    psFree(covariances);
+
     // Now make a new readout to go in the new cell
     pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
@@ -1339 +1390 @@
     newReadout->mask   = mosaicMask;
     newReadout->variance = mosaicVariances;
+    newReadout->covariance = mosaicCovariances;
     psFree(newReadout);                 // Drop reference
 

trunk/psModules/src/camera/pmReadoutFake.c

@@ -82 +82 @@
     }
     PS_ASSERT_PTR_NON_NULL(psf, false);
-    if (radius > 0 && isfinite(minFlux) && minFlux > 0.0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot define both minimum flux and fixed radius.");
-        return false;
-    }
 
     readout->image = psImageRecycle(readout->image, numCols, numRows, PS_TYPE_F32);
@@ -128 +124 @@
         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
         fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
-        float fakeRadius = radius > 0 ? radius :
-            PS_MAX(1.0, fakeModel->modelRadius(fakeModel->params, minFlux)); // Radius of fake source
+        float fakeRadius = 1.0;         // Radius of fake source
+        if (isfinite(minFlux)) {
+            fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
+        }
+        if (radius > 0) {
+            fakeRadius = PS_MAX(fakeRadius, radius);
+        }
 
         if (xOffset) {

trunk/psModules/src/camera/pmReadoutFake.h

@@ -13 +13 @@
 #include <pmSourceMasks.h>
 
+/// Generate a fake readout from vectors
+bool pmReadoutFakeFromVectors(pmReadout *readout, ///< Output readout
+                              int numCols, int numRows, ///< Dimension of image
+                              const psVector *x, const psVector *y, ///< Source coordinates
+                              const psVector *mag, ///< Source magnitudes
+                              const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL
+                              const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL
+                              const pmPSF *psf, ///< PSF for sources
+                              float minFlux, ///< Minimum flux to bother about; for setting source radius
+                              int radius, ///< Fixed radius for sources
+                              bool circularise, ///< Circularise PSF model?
+                              bool normalisePeak ///< Normalise the peak value?
+    );
+
 /// Generate a fake readout from an array of sources
-bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout, or NULL
+bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout
                               int numCols, int numRows, ///< Dimension of image
                               const psArray *sources, ///< Array of pmSource

trunk/psModules/src/objects/Makefile.am

@@ -53 +53 @@
         pmGrowthCurveGenerate.c \
         pmGrowthCurve.c \
-        pmSourceMatch.c
+        pmSourceMatch.c \
+        pmDetEff.c
 
 EXTRA_DIST = \
@@ -90 +91 @@
         pmTrend2D.h \
         pmGrowthCurve.h \
-        pmSourceMatch.h
+        pmSourceMatch.h \
+        pmDetEff.h
 
 CLEANFILES = *~

trunk/psModules/src/objects/pmSourceIO.c

@@ -42 +42 @@
 #include "pmSource.h"
 #include "pmModelClass.h"
+#include "pmDetEff.h"
 #include "pmSourceIO.h"
 
 #define BLANK_HEADERS "BLANK.HEADERS"   // Name of metadata in camera configuration containing header names
                                         // for putting values into a blank PHU
+
+// lookup the EXTNAME values used for table data and image header segments
+static bool sourceExtensions(psString *headname, // Extension name for header
+                             psString *dataname, // Extension name for data
+                             psString *deteffname, // Extension name for detection efficiency
+                             psString *xsrcname, // Extension name for extended sources
+                             psString *xfitname, // Extension name for extended fits
+                             const pmFPAfile *file, // File of interest
+                             const pmFPAview *view // View to level of interest
+                             )
+{
+    bool status;                        // Status of MD lookup
+
+    // Menu of EXTNAME rules
+    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
+    if (!menu) {
+        psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+        return false;
+    }
+
+    // EXTNAME for image header
+    if (headname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *headname = pmFPAfileNameFromRule(rule, file, view);
+    }
+
+    // EXTNAME for table data
+    if (dataname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *dataname = pmFPAfileNameFromRule(rule, file, view);
+    }
+
+    // EXTNAME for detection efficiency
+    if (deteffname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.DETEFF");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DETEFF in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *deteffname = pmFPAfileNameFromRule(rule, file, view);
+    }
+
+    // EXTNAME for extended source data table
+    if (xsrcname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *xsrcname = pmFPAfileNameFromRule (rule, file, view);
+    }
+
+    if (xfitname) {
+        // EXTNAME for extended source data table
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *xfitname = pmFPAfileNameFromRule (rule, file, view);
+    }
+
+    return true;
+}
+
 
 // translations between psphot object types and dophot object types
@@ -271 +345 @@
 
     char *exttype  = NULL;
-    char *dataname = NULL;
-    char *xsrcname = NULL;
-    char *xfitname = NULL;
-    char *headname = NULL;
 
     // if sources is NULL, write out an empty table
@@ -354 +424 @@
 
         // define the EXTNAME values for the different data segments:
-        {
-            // lookup the EXTNAME values used for table data and image header segments
-            char *rule = NULL;
-
-            // Menu of EXTNAME rules
-            psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
-            if (!menu) {
-                psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
-                return false;
-            }
-
-            // EXTNAME for image header
-            rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
-            if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
-                return false;
-            }
-            headname = pmFPAfileNameFromRule (rule, file, view);
-
-            // EXTNAME for table data
-            rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
-            if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
-                return false;
-            }
-            dataname = pmFPAfileNameFromRule (rule, file, view);
-
-            if (XSRC_OUTPUT) {
-              // EXTNAME for extended source data table
-              rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
-              if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config");
-                return false;
-              }
-              xsrcname = pmFPAfileNameFromRule (rule, file, view);
-            }
-            if (XFIT_OUTPUT) {
-              // EXTNAME for extended source data table
-              rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
-              if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config");
-                return false;
-              }
-              xfitname = pmFPAfileNameFromRule (rule, file, view);
-            }
+        psString headname = NULL;
+        psString dataname = NULL;
+        psString deteffname = NULL;
+        psString xsrcname = NULL;
+        psString xfitname = NULL;
+        if (!sourceExtensions(&headname, &dataname, &deteffname, XSRC_OUTPUT ? &xsrcname : NULL,
+                              XFIT_OUTPUT ? &xfitname : NULL, file, view)) {
+            return false;
         }
 
@@ -480 +513 @@
 
             // XXX these are case-sensitive since the EXTYPE is case-sensitive
-            status = false;
+            status = true;
             if (!strcmp (exttype, "SMPDATA")) {
-                status = pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_DEV_0")) {
-                status = pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_DEV_1")) {
-                status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_CAL_0")) {
-                status = pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_V1")) {
-                status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_V2")) {
-                status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
-            }
+                status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
+            }
+
+            if (deteffname) {
+                status &= pmReadoutWriteDetEff(file->fits, readout, outhead, deteffname);
+            }
+
             if (xsrcname) {
               if (!strcmp (exttype, "PS1_DEV_1")) {
-                  status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
+                  status &= pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
               }
               if (!strcmp (exttype, "PS1_CAL_0")) {
-                  status = pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
+                  status &= pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
               }
               if (!strcmp (exttype, "PS1_V1")) {
-                  status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
+                  status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
               }
               if (!strcmp (exttype, "PS1_V2")) {
-                  status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
+                  status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
               }
             }
             if (xfitname) {
               if (!strcmp (exttype, "PS1_DEV_1")) {
-                  status = pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname);
+                  status &= pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname);
               }
               if (!strcmp (exttype, "PS1_CAL_0")) {
-                  status = pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
+                  status &= pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
               }
               if (!strcmp (exttype, "PS1_V1")) {
-                  status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
+                  status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
               }
               if (!strcmp (exttype, "PS1_V2")) {
-                  status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
+                  status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
               }
             }
@@ -572 +610 @@
     // not needed if only one chip
     if (file->fpa->chips->n == 1) {
-        pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
-        return true;
+        pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
+        return true;
     }
 
@@ -885 +923 @@
         hdu = pmFPAviewThisHDU (view, file->fpa);
 
-        // lookup the EXTNAME values used for table data and image header segments
-        char *rule = NULL;
-        // Menu of EXTNAME rules
-        psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
-        if (!menu) {
-            psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
-            return false;
-        }
-        // EXTNAME for image header
-        rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
-        if (!rule) {
-            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
-            return false;
-        }
-        char *headname = pmFPAfileNameFromRule (rule, file, view);
-        // EXTNAME for table data
-        rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
-        if (!rule) {
-            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
-            return false;
-        }
-        char *dataname = pmFPAfileNameFromRule (rule, file, view);
+        // define the EXTNAME values for the different data segments:
+        psString headname = NULL;
+        psString dataname = NULL;
+        psString deteffname = NULL;
+        if (!sourceExtensions(&headname, &dataname, &deteffname, NULL, NULL, file, view)) {
+            return false;
+        }
 
         // advance to the IMAGE HEADER extension
@@ -958 +981 @@
                 sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
             }
+
+            if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {
+                psError(PS_ERR_IO, false, "Unable to read detection efficiency");
+                return false;
+            }
         }
 
@@ -1070 +1098 @@
 }
 
-    
+

trunk/psModules/src/psmodules.h

@@ -133 +133 @@
 #include <pmSourceVisual.h>
 #include <pmSourceMatch.h>
+#include <pmDetEff.h>
 
 // The following headers are from random locations, here because they cross bounds

trunk/psphot/src/Makefile.am

@@ -128 +128 @@
         psphotCheckStarDistribution.c  \
         psphotThreadTools.c            \
-        psphotAddNoise.c
+        psphotAddNoise.c               \
+        psphotEfficiency.c
 
 # dropped? psphotGrowthCurve.c

trunk/psphot/src/psphot.h

@@ -49 +49 @@
 bool            psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
 bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
-bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
+bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final);
 bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
 
 bool            psphotReplaceAllSources (psArray *sources, psMetadata *recipe);
-bool            psphotRemoveAllSources (psArray *sources, psMetadata *recipe);
+bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
 
 bool            psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
@@ -64 +64 @@
 bool            psphotGuessModel_Threaded (psThreadJob *job);
 
-bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf);
+bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf);
 bool            psphotMagnitudes_Threaded (psThreadJob *job);
 
@@ -76 +76 @@
 bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
+bool            psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources);
 
 // thread-related:
@@ -147 +148 @@
 pmPSF          *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe);
 bool            psphotSetHeaderNstars (psMetadata *recipe, psArray *sources);
-bool            psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe);
-bool            psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe);
-bool            psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add);
+bool            psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe);
+bool            psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe);
+bool            psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add);
 bool            psphotRadialPlot (int *kapa, const char *filename, pmSource *source);
 bool            psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe);

trunk/psphot/src/psphotAddNoise.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) {
+bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
   return psphotAddOrSubNoise (readout, sources, recipe, true);
 }
 
-bool psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) {
+bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
   return psphotAddOrSubNoise (readout, sources, recipe, false);
 }
 
-bool psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) {
+bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) {
 
     bool status = false;
@@ -42 +42 @@
     PS_ASSERT (status, false);
     if (isfinite(GAIN)) {
-        FACTOR /= GAIN;
+        FACTOR /= GAIN;
     }
 
@@ -70 +70 @@
         oldshape.sxy = PAR[PM_PAR_SXY];
 
-        // XXX can this be done more intelligently?
-        if (oldI0 == 0.0) continue;
-        if (!isfinite(oldI0)) continue;
+        // XXX can this be done more intelligently?
+        if (oldI0 == 0.0) continue;
+        if (!isfinite(oldI0)) continue;
 
         // increase size and height of source

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -12 +12 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
+bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) {
 
     bool status;

trunk/psphot/src/psphotMagnitudes.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
+bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf) {
 
     bool status = false;
@@ -64 +64 @@
 
             psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, psf);
+            psArrayAdd(job->args, 1, (pmPSF*)psf);    // Casting away const
             psArrayAdd(job->args, 1, binning);
             psArrayAdd(job->args, 1, backModel);

trunk/psphot/src/psphotReadout.c

@@ -61 +61 @@
     // display the backsub and backgnd images
     psphotVisualShowBackground (config, view, readout);
-    
+
     // run a single-model test if desired (exits from here if test is run)
     psphotModelTest (config, view, recipe);
@@ -224 +224 @@
     psphotMagnitudes(config, readout, view, sources, psf);
 
+    if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
+        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
+        psErrorClear();
+    }
+
     // replace failed sources?
     // psphotReplaceUnfitSources (sources);

trunk/psphot/src/psphotReadoutMinimal.c

@@ -97 +97 @@
     psphotMagnitudes(config, readout, view, sources, psf);
 
+    if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
+        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
+        psErrorClear();
+    }
+
     // drop the references to the image pixels held by each source
     psphotSourceFreePixels (sources);

trunk/psphot/src/psphotReplaceUnfit.c

@@ -47 +47 @@
 }
 
-bool psphotRemoveAllSources (psArray *sources, psMetadata *recipe) {
+bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) {
 
     bool status;

trunk/psphot/src/psphotSignificanceImage.c

@@ -22 +22 @@
     }
 
-    bool status_x, status_y;
-    float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
-    float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
-    if (status_x && status_y) {
-      // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?)
-      SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
-      NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
-      guess = false;
+    bool statusMajor, statusMinor;
+    float fwhmMajor = psMetadataLookupF32(&statusMajor, recipe, "FWHM_MAJ");
+    float fwhmMinor = psMetadataLookupF32(&statusMinor, recipe, "FWHM_MIN");
+    if (statusMajor && statusMinor) {
+        // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?)
+        if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) {
+            psWarning("fwhmMajor (%f) or fwhmMinor (%f) is bad!", fwhmMajor, fwhmMinor);
+        }
+        SIGMA_SMTH  = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrt(2.0*log(2.0)));
+        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+        guess = false;
     } else {
-      // if we do not know the FWHM, use the guess smoothing kernel supplied.
-      // it is a configuration error if these are not supplied
-      SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
-      PS_ASSERT (status, NULL);
-      NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
-      PS_ASSERT (status, NULL);
-      guess = true;
+        // if we do not know the FWHM, use the guess smoothing kernel supplied.
+        // it is a configuration error if these are not supplied
+        SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
+        PS_ASSERT (status, NULL);
+        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+        PS_ASSERT (status, NULL);
+        guess = true;
     }
     // record the actual smoothing sigma