Changeset 25383 for trunk/psModules

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)
• 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)

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