Changeset 25407

Timestamp:

Sep 15, 2009, 4:03:13 PM (17 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/eam_branches/20090715/psModules/src

Files:

• astrom/pmAstrometryWCS.c (modified) (6 diffs)
• astrom/pmAstrometryWCS.h (modified) (1 diff)
• camera/pmFPA.h (modified) (1 diff)
• camera/pmFPAMaskWeight.c (modified) (10 diffs)
• camera/pmFPAMaskWeight.h (modified) (6 diffs)
• camera/pmFPAMosaic.c (modified) (13 diffs)
• camera/pmReadoutFake.c (modified) (7 diffs, 1 prop)
• camera/pmReadoutFake.h (modified) (1 diff)
• concepts/pmConcepts.c (modified) (2 diffs)
• concepts/pmConceptsStandard.c (modified) (2 diffs)
• imcombine/pmPSFEnvelope.c (modified) (2 diffs)
• imcombine/pmStack.c (modified) (1 diff)
• imcombine/pmSubtraction.c (modified) (3 diffs)
• imcombine/pmSubtraction.h (modified) (2 diffs)
• imcombine/pmSubtractionAnalysis.c (modified) (9 diffs)
• imcombine/pmSubtractionAnalysis.h (modified) (1 diff)
• imcombine/pmSubtractionKernels.c (modified) (1 diff)
• imcombine/pmSubtractionKernels.h (modified) (1 diff)
• imcombine/pmSubtractionMatch.c (modified) (17 diffs)
• imcombine/pmSubtractionMatch.h (modified) (1 diff)
• imcombine/pmSubtractionStamps.c (modified) (1 diff)
• imcombine/pmSubtractionStamps.h (modified) (1 diff)
• objects/Makefile.am (modified) (2 diffs)
• objects/pmDetEff.c (copied) (copied from trunk/psModules/src/objects/pmDetEff.c )
• objects/pmDetEff.h (copied) (copied from trunk/psModules/src/objects/pmDetEff.h )
• objects/pmSourceIO.c (modified) (8 diffs)
• objects/pmSourceMatch.c (modified) (1 diff)
• psmodules.h (modified) (1 diff)

branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.c

@@ -289 +289 @@
     // test the CDELTi varient
     if (pcKeys) {
+        wcs->wcsCDkeys = 0;
         wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1");
         wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2");
@@ -334 +335 @@
     // test the CDi_j varient
     if (cdKeys) {
+        wcs->wcsCDkeys = 1;
+
         wcs->trans->x->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD1_1"); // == PC1_1
         wcs->trans->x->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD1_2"); // == PC1_2
@@ -375 +378 @@
     // XXX make it optional to write out CDi_j terms, or other versions
     // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity
-    double cdelt1 = wcs->cdelt1;
-    double cdelt2 = wcs->cdelt2;
-    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
-    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
-
-    // test the PC00i00j varient:
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
-    psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
-
-    // Elixir-style polynomial terms
-    // XXX currently, Elixir/DVO cannot accept mixed orders
-    // XXX need to respect the masks
-    // XXX is wcs->cdelt1,2 always consistent?
-    int fitOrder = wcs->trans->x->nX;
-    if (fitOrder > 1) {
+    if (!(wcs->wcsCDkeys)) {
+
+      double cdelt1 = wcs->cdelt1;
+      double cdelt2 = wcs->cdelt2;
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
+      
+      // test the PC00i00j varient:
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
+      psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
+      
+      // Elixir-style polynomial terms
+      // XXX currently, Elixir/DVO cannot accept mixed orders
+      // XXX need to respect the masks
+      // XXX is wcs->cdelt1,2 always consistent?
+      int fitOrder = wcs->trans->x->nX;
+      if (fitOrder > 1) {
         for (int i = 0; i <= fitOrder; i++) {
-            for (int j = 0; j <= fitOrder; j++) {
-                if (i + j < 2)
-                    continue;
-                if (i + j > fitOrder)
-                    continue;
-                sprintf (name, "PCA1X%1dY%1d", i, j);
-                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
-                sprintf (name, "PCA2X%1dY%1d", i, j);
-                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
-            }
+          for (int j = 0; j <= fitOrder; j++) {
+            if (i + j < 2)
+              continue;
+            if (i + j > fitOrder)
+              continue;
+            sprintf (name, "PCA1X%1dY%1d", i, j);
+            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
+            sprintf (name, "PCA2X%1dY%1d", i, j);
+            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
+          }
         }
         psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
-    }
-
-    // remove any existing 'CDi_j style' wcs keywords
-    if (psMetadataLookup(header, "CD1_1")) {
+      }
+      
+      // remove any existing 'CDi_j style' wcs keywords
+      if (psMetadataLookup(header, "CD1_1")) {
         psMetadataRemoveKey(header, "CD1_1");
         psMetadataRemoveKey(header, "CD1_2");
         psMetadataRemoveKey(header, "CD2_1");
         psMetadataRemoveKey(header, "CD2_2");
+      }
+    }
+    if (wcs->wcsCDkeys) {
+
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_2", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1]);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_1", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0]);
+      psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_2", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1]);
+
+      if (psMetadataLookup(header, "PC001001")) {
+        psMetadataRemoveKey(header, "PC001001");
+        psMetadataRemoveKey(header, "PC001002");
+        psMetadataRemoveKey(header, "PC002001");
+        psMetadataRemoveKey(header, "PC002002");
+      }
     }
 
@@ -535 +555 @@
         fpa->toSky->R -= 2.0*M_PI;
 
+    fpa->wcsCDkeys = wcs->wcsCDkeys;
+
     psTrace ("psastro", 5, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
              chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
@@ -682 +704 @@
     wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
 
+    wcs->wcsCDkeys = fpa->wcsCDkeys;
     psFree (toTPA);
 
@@ -922 +945 @@
     wcs->trans = psPlaneTransformAlloc (nXorder, nYorder);
     wcs->toSky = NULL;
+    wcs->wcsCDkeys = 0;
 
     memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE);

branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.h

@@ -23 +23 @@
     double crpix1, crpix2;
     double cdelt1, cdelt2;
+    bool wcsCDkeys;
     psProjection *toSky;
     psPlaneTransform *trans;

branches/eam_branches/20090715/psModules/src/camera/pmFPA.h

@@ -48 +48 @@
     psPlaneTransform *toTPA;  ///< Transformation from focal plane to tangent plane, or NULL
     psProjection *toSky;         ///< Projection from tangent plane to sky, or NULL
+    bool wcsCDkeys;
     // Information
     psMetadata *concepts;               ///< FPA-level concepts

branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.c

@@ -111 +111 @@
         // psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n");
         // return false;
-        psWarning("CELL.SATURATION is not set --- completely masking cell.\n");
-        saturation = NAN;
+        psWarning("CELL.SATURATION is not set --- completely masking cell.\n");
+        saturation = NAN;
     }
     float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
@@ -118 +118 @@
         // psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n");
         // return false;
-        psWarning("CELL.BAD is not set --- completely masking cell.\n");
-        bad = NAN;
+        psWarning("CELL.BAD is not set --- completely masking cell.\n");
+        bad = NAN;
     }
     psTrace("psModules.camera", 5, "Saturation: %f, bad: %f\n", saturation, bad);
 
-    // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN; 
+    // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN;
     // in this case, we have to set the mask as well
     float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
@@ -140 +140 @@
     // completely mask if SATURATION or BAD are invalid
     if (isnan(saturation) || isnan(bad) || isnan(gain) || isnan(readnoise)) {
-        psImageInit(mask, badMask);
-        return true;
+        psImageInit(mask, badMask);
+        return true;
     }
 
@@ -230 +230 @@
         // return false;
         psWarning("CELL.GAIN is not set --- setting variance to NAN\n");
-        gain = NAN;
+        gain = NAN;
     }
     float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
@@ -237 +237 @@
         // return false;
         psWarning("CELL.READNOISE is not set --- setting variance to NAN\n");
-        readnoise = NAN;
+        readnoise = NAN;
     }
     // if we have a non-NAN readnoise, then we need to ensure it has been updated (not necessary if NAN)
@@ -248 +248 @@
     if (isnan(gain) || isnan(readnoise)) {
         if (!readout->variance) {
-            // generate the image if needed
+            // generate the image if needed
             readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
         }
-        // XXX need to set the mask, if defined
+        // XXX need to set the mask, if defined
         psImageInit(readout->variance, NAN);
-        return true;
+        return true;
     }
 
@@ -262 +262 @@
 
         // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1.
-        // XXX this calculation is wrong: limit is 1 e-, but this is in DN
+        // XXX this calculation is wrong: limit is 1 e-, but this is in DN
         readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs");
         readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max",
@@ -276 +276 @@
     // apply a supplied readnoise map (NOTE: in DN, not electrons):
     if (noiseMap) {
-        psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
-        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
-        psFree (rdVar);
+        psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
+        psFree (rdVar);
     } else {
-        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
     }
 
@@ -362 +362 @@
 
 
-bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
-                                 psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
+bool pmReadoutVarianceRenormalise(const pmReadout *readout, psImageMaskType maskVal,
+                                  int sample, float minValid, float maxValid)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
@@ -371 +371 @@
     psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts
 
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    psStats *varianceStats = psStatsAlloc(meanStat);// Statistics for mean
-    if (!psImageBackground(varianceStats, NULL, variance, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for image");
-        psFree(varianceStats);
-        psFree(rng);
-        return false;
-    }
-    float meanVariance = varianceStats->robustMedian; // Mean variance
-    psFree(varianceStats);
-
-    psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean
-    if (!psImageBackground(imageStats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure stdev of image");
-        psFree(imageStats);
-        psFree(rng);
-        return false;
-    }
-    float stdevImage = imageStats->robustStdev; // Standard deviation of image
-    psFree(imageStats);
-    psFree(rng);
-
-    float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", correction);
-    psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32));
-
-    return true;
-}
-
-
-bool pmReadoutVarianceRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,
-                               psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
-{
-    PM_ASSERT_READOUT_NON_NULL(readout, false);
-    PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_VARIANCE(readout, false);
-
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
-
-    // Measure background
-    psStats *bgStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);// Statistics for background
-    if (!psImageBackground(bgStats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
-        psFree(bgStats);
-        psFree(rng);
-        return false;
-    }
-    float bgMean = bgStats->robustMedian; // Background level
-    float bgNoise = bgStats->robustStdev; // Background standard deviation
-    psFree(bgStats);
-    psTrace("psModules.camera", 5, "Background is %f +/- %f\n", bgMean, bgNoise);
-
-
-    // Construct kernels for flux measurement
-    // We use N(0,width) and N(0,width/sqrt(2)) kernels, following psphotSignificanceImage.
-    float sigFactor = 4.0 * M_PI * PS_SQR(width); // Factor for conversion from im/wt ratio to significance
-    int size = RENORM_NUM_SIGMA, fullSize = 2 * size + 1; // Half-size and full size of Gaussian
-    psVector *gauss = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian for weighting
-    psVector *gauss2 = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian squared
-    for (int i = 0, x = -size; i < fullSize; i++, x++) {
-        gauss->data.F32[i] = expf(-0.5 * PS_SQR(x) / PS_SQR(width));
-        gauss2->data.F32[i] = expf(-PS_SQR(x) / PS_SQR(width));
-    }
-
-    // Size of image
-    int numCols = image->numCols, numRows = image->numRows; // Size of images
-    int xSize = numCols - fullSize, ySize = numRows - fullSize; // Size of consideration
-    int xOffset = size, yOffset = size;       // Offset to region of consideration
-
-    // Measure fluxes
-    float peakFlux = RENORM_PEAK * bgNoise;     // Peak flux for fake sources
-    psVector *noise = psVectorAlloc(num, PS_TYPE_F32); // Measurements of the noise
-    psVector *source = psVectorAlloc(num, PS_TYPE_F32); // Measurements of fake sources
-    psVector *guess = psVectorAlloc(num, PS_TYPE_F32); // Guess at significance
-    psVector *photMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for fluxes
-    for (int i = 0; i < num; i++) {
-        // Coordinates of interest
-        int xPix = psRandomUniform(rng) * xSize + xOffset + 0.5;
-        int yPix = psRandomUniform(rng) * ySize + yOffset + 0.5;
-        psAssert(xPix - size >= 0 && xPix + size < numCols &&
-                 yPix - size >= 0 && yPix + size < numRows,
-                 "Bad pixel position: %d,%d", xPix, yPix);
-
-        // Weighted aperture photometry
-        // This has the same effect as smoothing the image by the window function
-        float sumNoise = 0.0;       // Sum for noise measurement
-        float sumSource = 0.0;      // Sum for source measurement
-        float sumVariance = 0.0;      // Sum for variance measurement
-        float sumGauss = 0.0, sumGauss2 = 0.0; // Sums of Gaussian kernels
-        for (int v = 0, y = yPix - size; v < fullSize; v++, y++) {
-            float xSumNoise = 0.0;  // Sum for noise measurement in x
-            float xSumSource = 0.0; // Sum for source measurement in x
-            float xSumVariance = 0.0; // Sum for variance measurement in x
-            float xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x
-            float yGauss = gauss->data.F32[v]; // Value of Gaussian in y
-            for (int u = 0, x = xPix - size; u < fullSize; u++, x++) {
-                if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int numPix = numCols * numRows;                         // Number of pixels
+    int num = PS_MAX(sample, numPix);                       // Number we care about
+    psVector *signoise = psVectorAllocEmpty(num, PS_TYPE_F32);   // Signal-to-noise values
+
+    if (num >= numPix) {
+        // We have an image smaller than Nsubset, so just loop over the image pixels
+        int index = 0;                  // Index for vector
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                    !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) {
                     continue;
                 }
-                float value = image->data.F32[y][x] - bgMean; // Value of image
-                float xGauss = gauss->data.F32[u]; // Value of Gaussian in x
-                float xGauss2 = gauss2->data.F32[u]; // Value of Gaussian^2 in x
-                xSumNoise += value * xGauss;
-                xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss;
-                xSumVariance += variance->data.F32[y][x] * xGauss2;
-                xSumGauss += xGauss;
-                xSumGauss2 += xGauss2;
-            }
-            float yGauss2 = gauss2->data.F32[v]; // Value of Gaussian^2 in y
-            sumNoise += xSumNoise * yGauss;
-            sumSource += xSumSource * yGauss;
-            sumVariance += xSumVariance * yGauss2;
-            sumGauss += xSumGauss * yGauss;
-            sumGauss2 += xSumGauss2 * yGauss2;
-        }
-
-        photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) &&
-                                                isfinite(sumVariance) && sumGauss > 0 && sumGauss2 > 0) ?
-                                               0 : 0xFF);
-
-        float smoothImageNoise = sumNoise / sumGauss; // Value of smoothed image pixel for noise
-        float smoothImageSource = sumSource / sumGauss; // Value of smoothed image pixel for source
-        float smoothVariance = sumVariance / sumGauss2; // Value of smoothed variance pixel
-
-        noise->data.F32[i] = smoothImageNoise;
-        source->data.F32[i] = smoothImageSource;
-        guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothVariance : 0.0;
-        psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n",
-                i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance);
-    }
-    psFree(gauss);
-    psFree(gauss2);
-    psFree(rng);
-
-    // Standard deviation of fluxes gives us the real significance
-    psStats *stdevStats = psStatsAlloc(stdevStat); // Statistics
-    if (!psVectorStats(stdevStats, noise, NULL, photMask, 0xFF)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure standard deviation of fluxes");
-        psFree(stdevStats);
-        psFree(noise);
-        psFree(source);
-        psFree(guess);
-        psFree(photMask);
+
+                signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]);
+                index++;
+            }
+        }
+        signoise->n = index;
+    } else {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+        int index = 0;                  // Index for vector
+        for (long i = 0; i < num; i++) {
+            // Pixel coordinates
+            int pixel = numPix * psRandomUniform(rng);
+            int x = pixel % numCols;
+            int y = pixel / numCols;
+
+            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) {
+                continue;
+            }
+
+            signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]);
+            index++;
+        }
+        signoise->n = index;
+        psFree(rng);
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics
+
+    if (!psVectorStats(stats, signoise, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics on S/N image");
+        psFree(signoise);
         return false;
     }
-    float stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of fluxes
-    psFree(stdevStats);
-    psFree(noise);
-    psTrace("psModules.camera", 5, "Standard deviation of fluxes is %f\n", stdev);
-
-    // Ratio of measured significance to guessed significance
-    psVector *ratio = psVectorAlloc(num, PS_TYPE_F32); // Ratio of measured to guess
-    for (int i = 0; i < num; i++) {
-        float measuredSig = PS_SQR(source->data.F32[i] / stdev); // Measured significance
-        ratio->data.F32[i] = measuredSig / guess->data.F32[i];
-        if (guess->data.F32[i] <= 0.0 || source->data.F32[i] <= 0.0 || !isfinite(ratio->data.F32[i])) {
-            photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF;
-        }
-        psTrace("psModules.camera", 9, "Ratio %d: %f, %f, %f\n",
-                i, guess->data.F32[i], measuredSig, ratio->data.F32[i]);
-    }
-    psFree(source);
-    psFree(guess);
-
-    psStats *meanStats = psStatsAlloc(meanStat | stdevStat); // Statistics
-    if (!psVectorStats(meanStats, ratio, NULL, photMask, 0xFF)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean ratio");
-        psFree(meanStats);
-        psFree(ratio);
-        psFree(photMask);
-        return false;
-    }
-    float meanRatio = psStatsGetValue(meanStats, meanStat); // Mean ratio
-    psTrace("psModules.camera", 5, "Mean significance ratio is %f +/- %f\n",
-            meanRatio, psStatsGetValue(meanStats, stdevStat));
-    psFree(meanStats);
-    psFree(ratio);
-    psFree(photMask);
-
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", meanRatio);
-    psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
+    psFree(signoise);
+
+    float covar = sqrtf(psImageCovarianceFactor(readout->covariance)); // Covariance factor
+    float correction = stats->robustStdev / covar; // Correction factor
+    psFree(stats);
+    psLogMsg("psModules.camera", PS_LOG_DETAIL, "Variance renormalisation factor is %f", correction);
+
+    // Check valid range of correction factor
+    if ((isfinite(minValid) && correction < minValid) || (isfinite(maxValid) && correction > maxValid)) {
+        psWarning("Variance renormalisation is outside valid range: %f vs %f:%f --- no correction made",
+                  correction, minValid, maxValid);
+        return true;
+    }
+
+    psBinaryOp(variance, variance, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
+
+    pmHDU *hdu = pmHDUFromReadout(readout); // HDU for readout
+    if (hdu)  {
+        psString history = NULL;
+        psStringAppend(&history, "Rescaled variance by %6.4f (stdev by %6.4f)",
+                       PS_SQR(correction), correction);
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
+        psFree(history);
+    }
 
     return true;
 }
 
-
-bool pmReadoutVarianceRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,
-                           psStatsOptions stdevStat, int width, psRandom *rng)
-{
-    PM_ASSERT_READOUT_NON_NULL(readout, false);
-    PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_VARIANCE(readout, false);
-    PS_ASSERT_INT_POSITIVE(width, false);
-
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
-    int numCols = image->numCols, numRows = image->numRows; // Size of images
-    int xNum = numCols / width + 1, yNum = numRows / width + 1; // Number of renormalisation regions
-    float xSize = numCols / (float)xNum, ySize = numRows / (float)yNum; // Size of renormalisation regions
-
-    psStats *meanStats = psStatsAlloc(meanStat), *stdevStats = psStatsAlloc(stdevStat); // Statistics
-    psVector *buffer = NULL;
-
-    for (int j = 0; j < yNum; j++) {
-        // Bounds in y
-        int yMin = j * ySize;
-        int yMax = (j + 1) * ySize;
-        for (int i = 0; i < xNum; i++) {
-            // Bounds in x
-            int xMin = i * xSize;
-            int xMax = (i + 1) * xSize;
-
-            psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
-            psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels
-            psImage *subVariance = psImageSubset(variance, region); // Sub image of the variance pixels
-            psImage *subMask = mask ? psImageSubset(mask, region) : NULL; // Sub-image of the mask pixels
-
-            if (!psImageBackground(stdevStats, &buffer, subImage, subMask, maskVal, rng) ||
-                !psImageBackground(meanStats, &buffer, subVariance, subMask, maskVal, rng)) {
-                // Nothing we can do about it, but don't want to keel over and die, so do our best to flag it.
-                psString regionStr = psRegionToString(region); // String with region
-                psWarning("Unable to measure statistics over %s", regionStr);
-                psFree(regionStr);
-                psErrorClear();
-                psImageInit(subVariance, NAN);
-                if (subMask) {
-                    psImageInit(subMask, maskVal);
-                }
-            } else {
-                double meanVar = psStatsGetValue(meanStats, meanStat); // Mean of variance map
-                double stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of image
-                psTrace("psModules.camera", 3,
-                        "Region [%d:%d,%d:%d] has variance %lf, but mean of variance map is %lf\n",
-                        xMin, xMax, yMin, yMax, PS_SQR(stdev), meanVar);
-                psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
-            }
-
-            psFree(subImage);
-            psFree(subVariance);
-            psFree(subMask);
-        }
-    }
-    psFree(meanStats);
-    psFree(stdevStats);
-    psFree(rng);
-    psFree(buffer);
-
-    return true;
-}
 
 

branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.h

@@ -15 +15 @@
 /// @addtogroup Camera Camera Layout
 /// @{
-
-#if 0
-/// Pixel mask values
-typedef enum {
-    PM_MASK_CLEAR    = 0x00,            ///< The pixel is not masked
-    PM_MASK_BLANK    = 0x01,            ///< The pixel is blank or has no (valid) data
-    PM_MASK_FLAT     = 0x02,            ///< The pixel is non-positive in the flat-field
-    PM_MASK_DETECTOR = 0x02,            ///< The detector pixel is bad (e.g., bad column, charge trap)
-    PM_MASK_SAT      = 0x04,            ///< The pixel is saturated in the image of interest
-    PM_MASK_BAD      = 0x04,            ///< The pixel is low in the image of interest
-    PM_MASK_RANGE    = 0x04,            ///< The pixel is out of range in the image of interest
-    PM_MASK_CR       = 0x08,            ///< The pixel is probably a CR
-    PM_MASK_SPARE1   = 0x10,            ///< Spare mask value
-    PM_MASK_SPARE2   = 0x20,            ///< Spare mask value
-    PM_MASK_SUSPECT  = 0x40,            ///< The pixel is suspected of being bad, but may not be
-    PM_MASK_MARK     = 0x80,            ///< The pixel is marked as temporarily ignored
-} pmMaskValue;
-#define PM_MASK_MARK 0x80            ///< The pixel is marked as temporarily ignored
-#define PM_MASK_SAT  0x04            ///< The pixel is saturated in the image of interest
-#endif
 
 /// Set a temporary readout mask using CELL.SATURATION and CELL.BAD
@@ -54 +34 @@
 /// can't be generated.
 bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance
-                          const psImage *noiseMap, ///< 2D image of the read noise in DN
+                          const psImage *noiseMap, ///< 2D image of the read noise in DN
                           bool poisson    ///< Include poisson variance (in addition to read noise)?
     );
@@ -73 +53 @@
 /// with HDU entry).  This is intended for most operations.
 bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance
-                          const psImage *noiseMap, ///< 2D image of the read noise in DN
+                          const psImage *noiseMap, ///< 2D image of the read noise in DN
                                bool poisson    ///< Include poisson variance (in addition to read noise)?
     );
@@ -83 +63 @@
                                    psImageMaskType sat, ///< Mask value to give saturated pixels
                                    psImageMaskType bad, ///< Mask value to give bad (low) pixels
-                                   const psImage *noiseMap, ///< 2D image of the read noise in DN
+                                   const psImage *noiseMap, ///< 2D image of the read noise in DN
                                    bool poisson ///< Include poisson variance (in addition to read noise)?
     );
@@ -93 +73 @@
                                 psImageMaskType sat, ///< Mask value to give saturated pixels
                                 psImageMaskType bad, ///< Mask value to give bad (low) pixels
-                                const psImage *noiseMap, ///< 2D image of the read noise in DN
+                                const psImage *noiseMap, ///< 2D image of the read noise in DN
                                 bool poisson ///< Include poisson variance (in addition to read noise)?
     );
@@ -100 +80 @@
 ///
 /// The variance map is adjusted so that the mean matches the actual pixel variance in the image
-bool pmReadoutVarianceRenormPixels(
-    const pmReadout *readout,           ///< Readout to normalise
-    psImageMaskType maskVal,                 ///< Value to mask
-    psStatsOptions meanStat,            ///< Statistic to measure the mean (of the variance map)
-    psStatsOptions stdevStat,           ///< Statistic to measure the stdev (of the image)
-    psRandom *rng                       ///< Random number generator
-    );
-
-/// Renormalise the variance map to match the actual photometry variance
-///
-/// The variance map is adjusted so that the actual significance of fake sources matches the
-/// guestimated significance
-bool pmReadoutVarianceRenormPhot(
+bool pmReadoutVarianceRenormalise(
     const pmReadout *readout,           ///< Readout to normalise
     psImageMaskType maskVal,            ///< Value to mask
-    int num,                            ///< Number of instances to measure over the image
-    float width,                        ///< Photometry width
-    psStatsOptions meanStat,            ///< Statistic to measure the mean
-    psStatsOptions stdevStat,           ///< Statistic to measure the stdev
-    psRandom *rng                       ///< Random number generator
-    );
-
-/// Renormalise the variance map to match the actual variance
-///
-/// The variance in the image is measured in patches, and the variance map is adjusted so that the mean for
-/// that patch corresponds.
-bool pmReadoutVarianceRenorm(const pmReadout *readout, // Readout to normalise
-                             psImageMaskType maskVal, // Value to mask
-                             psStatsOptions meanStat, // Statistic to measure the mean (of the variance map)
-                             psStatsOptions stdevStat, // Statistic to measure the stdev (of the image)
-                             int width,   // Width of patch (pixels)
-                             psRandom *rng // Random number generator (for sub-sampling images)
+    int sample,                         ///< Sample size
+    float minValid,                     ///< Minimum valid renormalisation, or NAN
+    float maxValid                      ///< Maximum valid renormalisation, or NAN
     );
 

branches/eam_branches/20090715/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
 

branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.c

@@ -28 +28 @@
 #define MODEL_TYPE "PS_MODEL_RGAUSS"    // Type of model to use
 #define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
-#define SOURCE_MASK (PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
 #define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
                     PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
@@ -50 +49 @@
 }
 
-bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
-                              const psVector *xOffset, const psVector *yOffset, const pmPSF *psf,
-                              float minFlux, int radius, bool circularise, bool normalisePeak)
+
+bool pmReadoutFakeFromVectors(pmReadout *readout, int numCols, int numRows,
+                              const psVector *x, const psVector *y, const psVector *mag,
+                              const psVector *xOffset, const psVector *yOffset,
+                              const pmPSF *psf, float minFlux, int radius,
+                              bool circularise, bool normalisePeak)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
     PS_ASSERT_INT_LARGER_THAN(numCols, 0, false);
     PS_ASSERT_INT_LARGER_THAN(numRows, 0, false);
-    PS_ASSERT_ARRAY_NON_NULL(sources, false);
-
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(mag, false);
+    PS_ASSERT_VECTOR_TYPE(mag, PS_TYPE_F32, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(mag, x, false);
+    long numSources = x->n;              // Number of sources
     if (xOffset || yOffset) {
         PS_ASSERT_VECTOR_NON_NULL(xOffset, false);
@@ -64 +73 @@
         PS_ASSERT_VECTORS_SIZE_EQUAL(xOffset, yOffset, false);
         PS_ASSERT_VECTOR_TYPE(xOffset, PS_TYPE_S32, false);
-        PS_ASSERT_VECTOR_TYPE_EQUAL(xOffset, yOffset, false);
-        if (xOffset->n != sources->n) {
+        PS_ASSERT_VECTOR_TYPE(yOffset, PS_TYPE_S32, false);
+        if (xOffset->n != numSources) {
             psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                     "Number of offset vectors (%ld) and sources (%ld) doesn't match",
-                    xOffset->n, sources->n);
+                    xOffset->n, numSources);
             return false;
         }
     }
     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);
     psImageInit(readout->image, 0);
 
-    int numSources = sources->n;          // Number of stars
-    for (int i = 0; i < numSources; i++) {
-        pmSource *source = sources->data[i]; // Source of interest
-        if (!source) {
-            continue;
-        }
-        if (source->mode & SOURCE_MASK) {
-            continue;
-        }
-        if (!isfinite(source->psfMag)) {
-            continue;
-        }
-        float x, y;                     // Coordinates of source
-        if (source->modelPSF) {
-            x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-            y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-        } else {
-            x = source->peak->xf;
-            y = source->peak->yf;
-        }
-
-        float flux = powf(10.0, -0.4 * source->psfMag); // Flux of source
+    for (long i = 0; i < numSources; i++) {
+        float flux = powf(10.0, -0.4 * mag->data.F32[i]); // Flux of source
+        float xSrc = x->data.F32[i], ySrc = y->data.F32[i]; // Coordinates of source
 
         if (normalisePeak) {
             // Normalise flux
-            pmModel *normModel = pmModelFromPSFforXY(psf, x, y, 1.0); // Model for normalisation
+            pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation
             if (!normModel || (normModel->flags & MODEL_MASK)) {
                 psFree(normModel);
@@ -121 +107 @@
         }
 
-        pmModel *fakeModel = pmModelFromPSFforXY(psf, x, y, flux);
+        pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux);
         if (!fakeModel || (fakeModel->flags & MODEL_MASK)) {
             psFree(fakeModel);
@@ -137 +123 @@
 
         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
-        fakeSource->peak = pmPeakAlloc(x, y, 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
+        fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
+        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) {
-            if (!pmSourceDefinePixels(fakeSource, readout, x + xOffset->data.S32[i],
-                                      y + yOffset->data.S32[i], fakeRadius)) {
+            if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[i],
+                                      ySrc + yOffset->data.S32[i], fakeRadius)) {
                 psErrorClear();
                 continue;
@@ -153 +144 @@
             }
         } else {
-            if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
+            if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) {
                 psErrorClear();
                 continue;
@@ -168 +159 @@
     return true;
 }
+
+
+bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
+                              pmSourceMode sourceMask, const psVector *xOffset, const psVector *yOffset,
+                              const pmPSF *psf, float minFlux, int radius,
+                              bool circularise, bool normalisePeak)
+{
+    PS_ASSERT_ARRAY_NON_NULL(sources, false);
+
+    int numSources = sources->n;          // Number of stars
+    psVector *x = psVectorAllocEmpty(numSources, PS_TYPE_F32);
+    psVector *y = psVectorAllocEmpty(numSources, PS_TYPE_F32);
+    psVector *mag = psVectorAllocEmpty(numSources, PS_TYPE_F32);
+
+    int numGood = 0;                    // Number of good sources
+    for (int i = 0; i < numSources; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+        if (source->mode & sourceMask) {
+            continue;
+        }
+        if (!isfinite(source->psfMag)) {
+            continue;
+        }
+        float xSrc, ySrc;                     // Coordinates of source
+        if (source->modelPSF) {
+            xSrc = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+            ySrc = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+        } else {
+            xSrc = source->peak->xf;
+            ySrc = source->peak->yf;
+        }
+
+        x->data.F32[numGood] = xSrc;
+        y->data.F32[numGood] = ySrc;
+        mag->data.F32[numGood] = source->psfMag;
+        numGood++;
+    }
+    x->n = numGood;
+    y->n = numGood;
+    mag->n = numGood;
+
+    bool status = pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf,
+                                           minFlux, radius, circularise, normalisePeak);
+    psFree(x);
+    psFree(y);
+    psFree(mag);
+
+    return status;
+}

branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.h

@@ -11 +11 @@
 #include <pmTrend2D.h>
 #include <pmPSF.h>
+#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
+                              pmSourceMode sourceMask, ///< Mask for sources
                               const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL
                               const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL

branches/eam_branches/20090715/psModules/src/concepts/pmConcepts.c

@@ -377 +377 @@
         conceptRegisterS32("CELL.XWINDOW", "Start of cell window (pixels)",p_pmConceptParse_Positions,p_pmConceptFormat_Positions, NULL, true, PM_FPA_LEVEL_CELL);
         conceptRegisterS32("CELL.YWINDOW", "Start of cell window (pixels)",p_pmConceptParse_Positions,p_pmConceptFormat_Positions, NULL, true, PM_FPA_LEVEL_CELL);
-        conceptRegisterF32("CELL.VARFACTOR", "Variance factor for conversion from large to small scales", NULL, NULL, NULL, true, PM_FPA_LEVEL_CELL);
 
         // CELL.TRIMSEC
@@ -485 +484 @@
         concept[length - 1] = '\0';
 
-        // special variants:
-        if (!strcmp(concept, "FPA.DATE")) {
-          psTime *fpaTime = psMetadataLookupPtr(NULL, fpa->concepts, "FPA.TIME");
-          if (!fpaTime) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Missing concept FPA.TIME needed for FPA.DATE");
-            psFree(string);
-            return NULL;
-          }
-          psString dateTimeString = psTimeToISO(fpaTime); // String representation
-          psList *dateTime = psStringSplit(dateTimeString, "T", true);
-          psFree(dateTimeString);
-          psString dateString = psMemIncrRefCounter(psListGet(dateTime, PS_LIST_HEAD)); // The date string
-          psFree (dateTime);
-
-          if (!psStringSubstitute(&string, dateString, "{FPA.DATE}")) {
-              psError(PS_ERR_UNKNOWN, false, "Unable to replace concept %s", concept);
-              psFree(string); 
-              psFree(dateString);
-              return NULL; 
-          }
-          psFree (dateString);
-          continue;
-        }
+        // special variants:
+        if (!strcmp(concept, "FPA.DATE")) {
+          psTime *fpaTime = psMetadataLookupPtr(NULL, fpa->concepts, "FPA.TIME");
+          if (!fpaTime) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Missing concept FPA.TIME needed for FPA.DATE");
+            psFree(string);
+            return NULL;
+          }
+          psString dateTimeString = psTimeToISO(fpaTime); // String representation
+          psList *dateTime = psStringSplit(dateTimeString, "T", true);
+          psFree(dateTimeString);
+          psString dateString = psMemIncrRefCounter(psListGet(dateTime, PS_LIST_HEAD)); // The date string
+          psFree (dateTime);
+
+          if (!psStringSubstitute(&string, dateString, "{FPA.DATE}")) {
+              psError(PS_ERR_UNKNOWN, false, "Unable to replace concept %s", concept);
+              psFree(string);
+              psFree(dateString);
+              return NULL;
+          }
+          psFree (dateString);
+          continue;
+        }
 
         psTrace("psModules.concepts", 7, "Interpolating concept %s", concept);

branches/eam_branches/20090715/psModules/src/concepts/pmConceptsStandard.c

@@ -429 +429 @@
         }
 
-        psString ra = psMetadataLookupStr(&mdok, formats, "FPA.RA"); // Format for RA
-        psString dec = psMetadataLookupStr(&mdok, formats, "FPA.DEC"); // Format for Dec
-        if (ra && strcasecmp(ra, "HOURS") == 0 && dec && strcasecmp(dec, "DEGREES") == 0) {
-            sexagesimal = true;
+        if (strcmp(concept->name, "FPA.RA") == 0 || strcmp(concept->name, "FPA.DEC") == 0) {
+            psString ra = psMetadataLookupStr(&mdok, formats, "FPA.RA"); // Format for RA
+            psString dec = psMetadataLookupStr(&mdok, formats, "FPA.DEC"); // Format for Dec
+            if (ra && strcasecmp(ra, "HOURS") == 0 && dec && strcasecmp(dec, "DEGREES") == 0) {
+                sexagesimal = true;
+            }
         }
     } else {
@@ -450 +452 @@
         small = 3600.0 * coords;
         small = (float)((int)(small * 1000.0)) / 1000.0;
-        if (negative) {
-            big *= -1;
-        }
         psString coordString = NULL;        // String with the coordinates in sexagesimal format
-        psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small);
+        psStringAppend(&coordString, "%s%02d:%02d:%06.3f",
+                       negative ? "-" : (strcmp(concept->name, "FPA.DEC") == 0 ? "+" : ""),
+                       big, medium, small);
         coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString);
         psFree(coordString);

branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c

@@ -124 +124 @@
         pmResiduals *resid = psf->residuals;// PSF residuals
         psf->residuals = NULL;
-        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, xOffset, yOffset, psf,
+        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf,
                                       NAN, radius, true, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout.");
@@ -298 +298 @@
         }
 
-        // measure the source moments: tophat windowing, no pixel S/N cutoff
+        // measure the source moments: tophat windowing, no pixel S/N cutoff
         if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) {
             // Can't do anything about it; limp along as best we can

branches/eam_branches/20090715/psModules/src/imcombine/pmStack.c

@@ -30 +30 @@
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
 #define PIXEL_MAP_BUFFER 2              // Number of entries to add to pixel map at a time
-#define ADD_VARIANCE                    // Allow additional variance (besides variance factor)?
+//#define ADD_VARIANCE                    // Allow additional variance (besides variance factor)?
 #define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.c

@@ -733 +733 @@
 
 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
-                              const psVector *deviations, psImage *subMask, float sigmaRej, int footprint)
+                              const psVector *deviations, psImage *subMask, float sigmaRej)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
@@ -821 +821 @@
     ds9num++;
 
+    int footprint = stamps->footprint;  // Half-size of stamp region of interest
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
@@ -956 +957 @@
     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
 
-    psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return
-    psVector *fakeSolution = psVectorAlloc(solution->n, PS_TYPE_F64); // Fake solution vector
+    psArray *images = psArrayAlloc(kernels->solution1->n - 1); // Images of each kernel to return
+    psVector *fakeSolution = psVectorAlloc(kernels->solution1->n, PS_TYPE_F64); // Fake solution vector
     psVectorInit(fakeSolution, 0.0);
 
-    for (int i = 0; i < solution->n - 1; i++) {
-        fakeSolution->data.F64[i] = solution->data.F64[i];
+    for (int i = 0; i < kernels->solution1->n - 1; i++) {
+        fakeSolution->data.F64[i] = kernels->solution1->data.F64[i];
         images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual);
         fakeSolution->data.F64[i] = 0.0;

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.h

@@ -68 +68 @@
                               const psVector *deviations, ///< Deviations for each stamp
                               psImage *subMask, ///< Subtraction mask
-                              float sigmaRej, ///< Number of RMS deviations above zero at which to reject
-                              int footprint ///< Half-size of stamp
+                              float sigmaRej ///< Number of RMS deviations above zero at which to reject
     );
 
@@ -94 +93 @@
 
 /// Generate images of the convolution kernel elements
-psArray *pmSubtractionKernelSolutions(const psVector *solution, ///< Solution vector
-                                      const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                      float x, float y ///< Normalised position [-1,1] for images
+psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                      float x, float y, ///< Normalised position [-1,1] for images
+                                      bool wantDual ///< Calculate for the dual kernel?
     );
+
 
 /// Execute a thread job to convolve a patch of the image

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.c

@@ -17 +17 @@
 
 
-bool pmSubtractionAnalysis(psMetadata *analysis, pmSubtractionKernels *kernels, psRegion *region,
+bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header,
+                           pmSubtractionKernels *kernels, psRegion *region,
                            int numCols, int numRows)
 {
-    if (analysis) {
-        PS_ASSERT_METADATA_NON_NULL(analysis, false);
-    }
+    PS_ASSERT_METADATA_NON_NULL(analysis, false);
+    PS_ASSERT_METADATA_NON_NULL(header, false);
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false);
@@ -38 +38 @@
                          PS_DATA_REGION | PS_META_DUPLICATE_OK,
                          "Region over which subtraction was performed", subRegion);
+
+        psString string = psRegionToString(*subRegion);
         psFree(subRegion);
+
+        psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK,
+                         "Region over which subtraction was performed", string);
+        psFree(string);
     }
 
@@ -45 +51 @@
                      PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
     psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
-                     PS_META_DUPLICATE_OK, "Subtraction kernels", kernels->mode);
+                     PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
+    psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
+                     PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
 
     // Realisations of kernel
@@ -113 +121 @@
     {
         psMetadata *header = psMetadataAlloc(); // Header
-        for (int i = 0; i < solution->n; i++) {
+        for (int i = 0; i < kernels->solution1->n; i++) {
             psString name = NULL;       // Header keyword
             psStringAppend(&name, "SOLN%04d", i);
-            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
+            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, kernels->solution1->data.F64[i]);
             psFree(name);
         }
-        psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
+        psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, false);
         psFits *kernelFile = psFitsOpen("kernels.fits", "w");
         (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
@@ -128 +136 @@
 #endif
 
-
-    // Set the variance factors
-    float vf1 = 1.0, vf2 = 1.0;         // Variance factors for each image
-    switch (kernels->mode) {
-      case PM_SUBTRACTION_MODE_1:
-        vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
-        break;
-      case PM_SUBTRACTION_MODE_2:
-        vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
-        break;
-      case PM_SUBTRACTION_MODE_DUAL:
-        vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
-        vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, true);
-        break;
-      default:
-        psAbort("Invalid subtraction mode: %x", kernels->mode);
-    }
-
-    // Weight by the area
-    if (region) {
-        float norm = (region->x1 - region->x0 + 1) * (region->y1 - region->y0 + 1) / (numCols * numRows);
-        vf1 *= norm;
-        vf2 *= norm;
-    }
-
-    // Update the variance factor
-#define UPDATE_VARFACTOR(VF, ANALYSIS) { \
-    psMetadataItem *vfItem = psMetadataLookup(analysis, ANALYSIS); \
-    if (vfItem) { \
-        psAssert(vfItem->type == PS_TYPE_F32, "Should be the type we said."); \
-        vfItem->data.F32 += VF; \
-    } else { \
-        psMetadataAddF32(analysis, PS_LIST_TAIL, ANALYSIS, 0, "Variance factor weighted by the area", VF); \
-    } \
-}
-
-    UPDATE_VARFACTOR(vf1, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
-    UPDATE_VARFACTOR(vf2, PM_SUBTRACTION_ANALYSIS_VARFACTOR_2);
 
     // Kernel shape
@@ -201 +171 @@
         psFree(image);
 
-        psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
-        if (item) {
-            item->data.F32 = PS_MAX(item->data.F32, max);
-        } else {
-            psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
-                             "Maximum deconvolution fraction", max);
+        {
+            psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
+            if (item) {
+                max = item->data.F32 = PS_MAX(item->data.F32, max);
+            } else {
+                psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
+                                 "Maximum deconvolution fraction", max);
+            }
+        }
+
+        {
+            psMetadataItem *item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
+            if (item) {
+                item->data.F32 = max;
+            } else {
+                psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
+                                 "Maximum deconvolution fraction", max);
+            }
         }
     }
@@ -254 +236 @@
         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY,
                          PS_META_DUPLICATE_OK, "Moment in yy", m02);
+
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM,
+                         PS_META_DUPLICATE_OK, "Normalisation", m00);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX,
+                         PS_META_DUPLICATE_OK, "Moment in x", m10);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MY,
+                         PS_META_DUPLICATE_OK, "Moment in y", m01);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXX,
+                         PS_META_DUPLICATE_OK, "Moment in xx", m20);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXY,
+                         PS_META_DUPLICATE_OK, "Moment in xy", m11);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY,
+                         PS_META_DUPLICATE_OK, "Moment in yy", m02);
     }
 
@@ -263 +258 @@
 
         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF,
+                         PS_META_DUPLICATE_OK, "Background difference", bg);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF,
                          PS_META_DUPLICATE_OK, "Background difference", bg);
         psFree(polyValues);
@@ -275 +272 @@
         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
                          kernels->rms);
+
+        psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps",
+                         kernels->numStamps);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation",
+                         kernels->mean);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
+                         kernels->rms);
     }
 

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.h

@@ -27 +27 @@
 bool pmSubtractionAnalysis(
     psMetadata *analysis,               ///< Metadata container for QA information
+    psMetadata *header,                 ///< Metadata container for QA information to put in header
     pmSubtractionKernels *kernels,      ///< Kernels
     psRegion *region,                   ///< Region for subtraction

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.c

@@ -765 +765 @@
     return PM_SUBTRACTION_KERNEL_NONE;
 }
+
+pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in)
+{
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(in, NULL);
+
+    pmSubtractionKernels *out = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
+    psMemSetDeallocator(out, (psFreeFunc)subtractionKernelsFree);
+
+    out->type = in->type;
+    out->description = psMemIncrRefCounter(in->description);
+    out->num = in->num;
+    out->u = psMemIncrRefCounter(in->u);
+    out->v = psMemIncrRefCounter(in->v);
+    out->widths = psMemIncrRefCounter(in->widths);
+    out->preCalc = psMemIncrRefCounter(in->preCalc);
+    out->penalty = in->penalty;
+    out->penalties = psMemIncrRefCounter(in->penalties);
+    out->uStop = psMemIncrRefCounter(in->uStop);
+    out->vStop = psMemIncrRefCounter(in->vStop);
+    out->size = in->size;
+    out->inner = in->inner;
+    out->spatialOrder = in->spatialOrder;
+    out->bgOrder = in->bgOrder;
+    out->mode = in->mode;
+    out->numCols = in->numCols;
+    out->numRows = in->numRows;
+    out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
+    out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
+
+    return out;
+}

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.h

@@ -203 +203 @@
     );
 
+/// Copy kernels
+///
+/// A deep copy is performed on the solution only; the other components are merely pointers.
+pmSubtractionKernels *pmSubtractionKernelsCopy(
+    const pmSubtractionKernels *in      // Kernels to copy
+    );
+
 
 #endif

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c

@@ -10 +10 @@
 #include "pmHDU.h"
 #include "pmFPA.h"
+#include "pmHDUUtils.h"
 #include "pmSubtractionParams.h"
 #include "pmSubtractionKernels.h"
@@ -58 +59 @@
 
 
-static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read
-                      const pmReadout *ro1, // Readout 1
-                      const pmReadout *ro2, // Readout 2
-                      const psImage *subMask, // Mask for subtraction, or NULL
-                      psImage *variance,  // Variance map
-                      const psRegion *region, // Region of interest, or NULL
-                      float thresh1,  // Threshold for stamp finding on readout 1
-                      float thresh2,  // Threshold for stamp finding on readout 2
-                      float stampSpacing, // Spacing between stamps
-                      int size,         // Kernel half-size
-                      int footprint,     // Convolution footprint for stamps
-                      pmSubtractionMode mode // Mode for subtraction
+static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read
+                                 const pmReadout *ro1, // Readout 1
+                                 const pmReadout *ro2, // Readout 2
+                                 const psImage *subMask, // Mask for subtraction, or NULL
+                                 psImage *variance,  // Variance map
+                                 const psRegion *region, // Region of interest, or NULL
+                                 float thresh1,  // Threshold for stamp finding on readout 1
+                                 float thresh2,  // Threshold for stamp finding on readout 2
+                                 float stampSpacing, // Spacing between stamps
+                                 int size,         // Kernel half-size
+                                 int footprint,     // Convolution footprint for stamps
+                                 pmSubtractionMode mode // Mode for subtraction
     )
 {
@@ -163 +164 @@
 }
 
+static void subtractionAnalysisUpdate(pmReadout *conv1, pmReadout *conv2, // Convolved images
+                                      const psMetadata *analysis, // Analysis metadata
+                                      const psMetadata *header // Header metadata
+    )
+{
+    if (conv1) {
+        conv1->analysis = psMetadataCopy(conv1->analysis, analysis);
+    }
+    if (conv2) {
+        conv2->analysis = psMetadataCopy(conv2->analysis, analysis);
+    }
+
+    if (conv1 && conv1->parent) {
+        pmHDU *hdu = pmHDUFromCell(conv1->parent);
+        if (hdu) {
+            hdu->header = psMetadataCopy(hdu->header, header);
+        }
+    }
+    if (conv2 && conv2->parent) {
+        pmHDU *hdu = pmHDUFromCell(conv2->parent);
+        if (hdu) {
+            hdu->header = psMetadataCopy(hdu->header, header);
+        }
+    }
+
+    return;
+}
 
 
@@ -253 +281 @@
 
     psMetadata *outAnalysis = psMetadataAlloc(); // Output analysis values
+    psMetadata *outHeader = psMetadataAlloc(); // Output header values
 
     psTrace("psModules.imcombine", 2, "Convolving...\n");
@@ -259 +288 @@
         psRegion *region = regions->data[i]; // Region of interest
 
-        if (!pmSubtractionAnalysis(outAnalysis, kernel, region, ro1->image->numCols, ro1->image->numRows)) {
+        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region,
+                                   ro1->image->numCols, ro1->image->numRows)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
             psFree(outAnalysis);
+            psFree(outHeader);
             psFree(subMask);
             psFree(kernels);
@@ -272 +303 @@
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
             psFree(outAnalysis);
+            psFree(outHeader);
             psFree(subMask);
             psFree(kernels);
@@ -283 +315 @@
     psFree(regions);
 
-    if (conv1) {
-        psMetadataCopy(conv1->analysis, outAnalysis);
-    }
-    if (conv2) {
-        psMetadataCopy(conv2->analysis, outAnalysis);
-    }
+    subtractionAnalysisUpdate(conv1, conv2, outAnalysis, outHeader);
     psFree(outAnalysis);
+    psFree(outHeader);
 
     return true;
@@ -369 +397 @@
     pmSubtractionKernels *kernels = NULL; // Kernel basis functions
     psMetadata *analysis = psMetadataAlloc(); // QA data
+    psMetadata *header = psMetadataAlloc(); // QA data for header
 
     int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
@@ -442 +471 @@
             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
             // doesn't matter.
-            if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
+            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
                            stampSpacing, size, footprint, subMode)) {
                 goto MATCH_ERROR;
             }
 
+
+            // Define kernel basis functions
+            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
+                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
+                                                          stamps, footprint, optThreshold, penalty, subMode);
+                if (!kernels) {
+                    psErrorClear();
+                    psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
+                }
+            }
+            if (kernels == NULL) {
+                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
+                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
+                                                       inner, binning, ringsOrder, penalty, subMode);
+            }
+
+            memCheck("kernels");
+
             if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
+#if 0
                 // Get backgrounds
                 psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background
@@ -469 +517 @@
 
                 pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use
+#endif
+                pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
                 switch (newMode) {
                   case PM_SUBTRACTION_MODE_1:
@@ -483 +533 @@
             }
 
-            // Define kernel basis functions
-            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
-                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                          stamps, footprint, optThreshold, penalty, subMode);
-                if (!kernels) {
-                    psErrorClear();
-                    psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
-                }
-            }
-            if (kernels == NULL) {
-                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
-                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
-                                                       inner, binning, ringsOrder, penalty, subMode);
-            }
-
-            memCheck("kernels");
-
             int numRejected = -1;       // Number of rejected stamps in each iteration
             for (int k = 0; k < iter && numRejected != 0; k++) {
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-                if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
-                               stampSpacing, size, footprint, subMode)) {
+                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
+                                          stampThresh1, stampThresh2, stampSpacing,
+                                          size, footprint, subMode)) {
                     goto MATCH_ERROR;
                 }
@@ -535 +569 @@
 
                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
-                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint);
+                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
                 if (numRejected < 0) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
@@ -557 +591 @@
                     goto MATCH_ERROR;
                 }
-                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint);
+                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
                 psFree(deviations);
             }
@@ -565 +599 @@
             memCheck("solution");
 
-            if (!pmSubtractionAnalysis(analysis, kernels, region, numCols, numRows)) {
+            if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
                 goto MATCH_ERROR;
@@ -601 +635 @@
     memCheck("convolution");
 
-    if (conv1) {
-        psMetadataCopy(conv1->analysis, analysis);
-    }
-    if (conv2) {
-        psMetadataCopy(conv2->analysis, analysis);
-    }
+    subtractionAnalysisUpdate(conv1, conv2, analysis, header);
     psFree(analysis);
+    psFree(header);
 
 #ifdef TESTING
@@ -629 +659 @@
 MATCH_ERROR:
     psFree(analysis);
+    psFree(header);
     psFree(region);
     psFree(regionString);
@@ -842 +873 @@
     return mode;
 }
+
+
+// Test a subtraction mode by performing a single iteration
+static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode
+                                pmSubtractionKernels *kernels, // Kernel description
+                                const char *description, // Description for trace
+                                psImage *subMask,  // Subtraction mask
+                                float rej               // Rejection threshold
+                                )
+{
+    assert(stamps);
+    assert(kernels);
+
+    psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description);
+    if (!pmSubtractionSolveEquation(kernels, stamps)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
+    psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
+    if (!deviations) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
+    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
+    if (numRejected < 0) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+        psFree(deviations);
+        return false;
+    }
+    psFree(deviations);
+
+    if (numRejected > 0) {
+        // Allow re-fit with reduced stamps set
+        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
+        if (!pmSubtractionSolveEquation(kernels, stamps)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
+        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
+        if (!deviations) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description);
+        long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
+        if (numRejected < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+            psFree(deviations);
+            return false;
+        }
+        psFree(deviations);
+    }
+
+    return true;
+}
+
+
+pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels,
+                                        const psImage *subMask, float rej)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(*kernels, PM_SUBTRACTION_MODE_ERR);
+
+    // Copies of the inputs
+    pmSubtractionStampList *stamps1 = pmSubtractionStampListCopy(*stamps);
+    pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(*kernels);
+    psImage *subMask1 = psImageCopy(NULL, subMask, subMask->type.type);
+    kernels1->mode = PM_SUBTRACTION_MODE_1;
+
+    if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 1");
+        psFree(stamps1);
+        psFree(kernels1);
+        psFree(subMask1);
+        return PM_SUBTRACTION_MODE_ERR;
+    }
+    psFree(subMask1);
+
+    // Copies of the inputs
+    pmSubtractionStampList *stamps2 = pmSubtractionStampListCopy(*stamps);
+    pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(*kernels);
+    psImage *subMask2 = psImageCopy(NULL, subMask, subMask->type.type);
+    kernels2->mode = PM_SUBTRACTION_MODE_2;
+
+    if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 2");
+        psFree(stamps2);
+        psFree(kernels2);
+        psFree(subMask2);
+        psFree(stamps1);
+        psFree(kernels1);
+        return PM_SUBTRACTION_MODE_ERR;
+    }
+    psFree(subMask2);
+
+
+    pmSubtractionStampList *bestStamps = NULL; // Best choice for stamps
+    pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels
+    psLogMsg("psModules.imcombine", PS_LOG_INFO,
+             "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n",
+             kernels1->mean, kernels1->rms, kernels1->numStamps,
+             kernels2->mean, kernels2->rms, kernels2->numStamps);
+
+    if (kernels1->mean < kernels2->mean) {
+        bestStamps = stamps1;
+        bestKernels = kernels1;
+    } else {
+        bestStamps = stamps2;
+        bestKernels = kernels2;
+    }
+
+    psFree(*stamps);
+    psFree(*kernels);
+    *stamps = psMemIncrRefCounter(bestStamps);
+    *kernels = psMemIncrRefCounter(bestKernels);
+
+    psFree(stamps1);
+    psFree(stamps2);
+    psFree(kernels1);
+    psFree(kernels2);
+
+    return bestKernels->mode;
+}
+

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.h

@@ -83 +83 @@
     );
 
+/// Determine best subtraction mode to use
+///
+/// Subtractions are attempted each way, and the mode with the lower residual is taken to be the best
+pmSubtractionMode pmSubtractionBestMode(
+    pmSubtractionStampList **stamps,    ///< Stamps to use for solution
+    pmSubtractionKernels **kernels,     ///< Kernels to use for solution
+    const psImage *subMask,             ///< Subtraction mask
+    float rej                           ///< Rejection threshold for stamps
+    );
 
 #endif

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.c

@@ -225 +225 @@
 }
 
+pmSubtractionStampList *pmSubtractionStampListCopy(const pmSubtractionStampList *in)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(in, NULL);
+
+    pmSubtractionStampList *out = psAlloc(sizeof(pmSubtractionStampList)); // Copied stamp list to return
+    psMemSetDeallocator(out, (psFreeFunc)subtractionStampListFree);
+
+    int num = out->num = in->num;       // Number of stamps
+    out->stamps = psArrayAlloc(num);
+    out->regions = psArrayAlloc(num);
+    out->x = NULL;
+    out->y = NULL;
+    out->flux = NULL;
+    out->footprint = in->footprint;
+
+    for (int i = 0; i < num; i++) {
+        psRegion *inRegion = in->regions->data[i]; // Input region
+        out->regions->data[i] = psRegionAlloc(inRegion->x0, inRegion->x1, inRegion->y0, inRegion->y1);
+
+        pmSubtractionStamp *inStamp = in->stamps->data[i]; // Input stamp
+        pmSubtractionStamp *outStamp = psAlloc(sizeof(pmSubtractionStamp));
+        psMemSetDeallocator(outStamp, (psFreeFunc)subtractionStampFree);
+        outStamp->x = inStamp->x;
+        outStamp->y = inStamp->y;
+        outStamp->flux = inStamp->flux;
+        outStamp->xNorm = inStamp->xNorm;
+        outStamp->yNorm = inStamp->yNorm;
+        outStamp->status = inStamp->status;
+
+        outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL;
+        outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL;
+        outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL;
+
+        if (inStamp->convolutions1) {
+            int size = inStamp->convolutions1->n; // Size of array
+            outStamp->convolutions1 = psArrayAlloc(size);
+            for (int j = 0; j < size; j++) {
+                psKernel *conv = inStamp->convolutions1->data[j]; // Convolution
+                outStamp->convolutions1->data[j] = conv ? psKernelCopy(conv) : NULL;
+            }
+        } else {
+            outStamp->convolutions1 = NULL;
+        }
+        if (inStamp->convolutions2) {
+            int size = inStamp->convolutions2->n; // Size of array
+            outStamp->convolutions2 = psArrayAlloc(size);
+            for (int j = 0; j < size; j++) {
+                psKernel *conv = inStamp->convolutions2->data[j]; // Convolution
+                outStamp->convolutions2->data[j] = conv ? psKernelCopy(conv) : NULL;
+            }
+        } else {
+            outStamp->convolutions2 = NULL;
+        }
+
+        outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL;
+        outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL;
+        outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL;
+        outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL;
+        outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL;
+
+        out->stamps->data[i] = outStamp;
+    }
+
+    return out;
+}
+
 pmSubtractionStamp *pmSubtractionStampAlloc(void)
 {

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.h

@@ -52 +52 @@
     } \
 }
+
+/// Copy a list of stamps
+///
+/// A deep copy is performed of the stamp list and the component stamps
+pmSubtractionStampList *pmSubtractionStampListCopy(
+    const pmSubtractionStampList *in    // Stamp list to copy
+    );
+
 
 /// A stamp for image subtraction

branches/eam_branches/20090715/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 = *~

branches/eam_branches/20090715/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 @@
 }
 
-    
+

branches/eam_branches/20090715/psModules/src/objects/pmSourceMatch.c

@@ -221 +221 @@
         } else {
             // Match with the master list
-            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
+            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, PS_TREE_EUCLIDEAN, xMaster, yMaster); // kd Tree
             long numMatch = 0;          // Number of matches
 

branches/eam_branches/20090715/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