Changeset 6957 for trunk/psModules/src/detrend/pmFringeStats.c

Timestamp:

Apr 21, 2006, 5:21:00 PM (20 years ago)

Author:

Paul Price

Message:

Fringe subtraction is working!

File:

• trunk/psModules/src/detrend/pmFringeStats.c (modified) (11 diffs)

trunk/psModules/src/detrend/pmFringeStats.c

@@ -3 +3 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-20 04:28:00 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-22 03:21:00 $
  *
  *  Copyright 2004 IfA
  */
 
+#include <stdio.h>
+#include <assert.h>
 #include "pslib.h"
 #include "pmFPA.h"
@@ -15 +17 @@
 #define MAX(x,y) ((x) > (y) ? (x) : (y))
 #define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+
+// Future optimisations for speed:
+//
+// 1. Clipping --- don't re-do the matrix setup again, but carry matrix and vector around, subtract
+// contributions from clipped data points.
+// 2. Faster psImageStats (use memcpy?)
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// pmFringeRegions
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 static void fringeRegionsFree(pmFringeRegions *fringe)
@@ -44 +59 @@
 }
 
-
-static void fringeStatsFree(pmFringeStats *stats)
-{
-    psFree(stats->regions);
-    psFree(stats->f);
-    psFree(stats->df);
-}
-
-pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions)
-{
-    pmFringeStats *stats = psAlloc(sizeof(pmFringeStats));
-    (void)psMemSetDeallocator(stats, (psFreeFunc)fringeStatsFree);
-
-    stats->regions = psMemIncrRefCounter(regions);
-    stats->f = psVectorAlloc(regions->x->n, PS_TYPE_F32);
-    stats->f->n = regions->x->n;
-    stats->df = psVectorAlloc(regions->x->n, PS_TYPE_F32);
-    stats->df->n = regions->x->n;
-
-    return stats;
-}
-
-
-static void fringeScaleFree(pmFringeScale *scale)
-{
-    psFree(scale->coeff);
-    psFree(scale->coeffErr);
-    return;
-}
-
-pmFringeScale *pmFringeScaleAlloc(int nFringeFrames)
-{
-    pmFringeScale *scale = psAlloc(sizeof(pmFringeScale));
-    (void)psMemSetDeallocator(scale, (psFreeFunc)fringeScaleFree);
-
-    scale->nFringeFrames = nFringeFrames;
-    scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
-    scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
-    scale->coeff->n = nFringeFrames + 1;
-    scale->coeff->n = nFringeFrames + 1;
-
-    return scale;
-}
-
-bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, psImage *image)
+bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, const psImage *image)
 {
 
@@ -100 +71 @@
     fringe->mask = psVectorRecycle(fringe->mask, fringe->nRequested, PS_TYPE_U8);
     fringe->x->n = fringe->y->n = fringe->mask->n = fringe->nRequested;
+    psVectorInit(fringe->mask, 0);
 
     int nX = image->numCols;
@@ -120 +92 @@
 }
 
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// pmFringeStats
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+static void fringeStatsFree(pmFringeStats *stats)
+{
+    psFree(stats->regions);
+    psFree(stats->f);
+    psFree(stats->df);
+}
+
+pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions)
+{
+    pmFringeStats *stats = psAlloc(sizeof(pmFringeStats));
+    (void)psMemSetDeallocator(stats, (psFreeFunc)fringeStatsFree);
+
+    int numRegions = regions->x->n;     // Number of regions
+    stats->regions = psMemIncrRefCounter(regions);
+    stats->f = psVectorAlloc(numRegions, PS_TYPE_F32);
+    stats->df = psVectorAlloc(numRegions, PS_TYPE_F32);
+    stats->f->n = stats->df->n = numRegions;
+
+    return stats;
+}
+
 pmFringeStats *pmFringeStatsMeasure(pmFringeRegions *fringe, pmReadout *readout, psMaskType maskVal)
 {
@@ -143 +141 @@
     psImage *mask  = readout->mask;
 
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Median statistics only
+    psStats *medianSd = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); // Median and SD
+
+    // Measure the sky over the image
+    psImage *sky = psImageAlloc(fringe->nX, fringe->nY, PS_TYPE_F32);
+    for (int i = 0; i < fringe->nY; i++) {
+        int y0 = (float)i * (float)image->numRows / (float)fringe->nY;
+        int y1 = (float)(i + 1) * (float)image->numRows / (float)fringe->nY;
+        for (int j = 0; j < fringe->nX; j++) {
+            int x0 = (float)j * (float)image->numCols / (float)fringe->nX;
+            int x1 = (float)(j + 1) * (float)image->numCols / (float)fringe->nX;
+            psRegion region = psRegionSet(x0, x1, y0, y1);
+            psImage *subImage = psImageSubset(image, region); // Subimage of the sky region
+            psImage *subMask = psImageSubset(mask, region); // Subimage of the sky region
+            psImageStats(median, subImage, subMask, maskVal);
+            sky->data.F32[i][j] = median->sampleMedian;
+        }
+    }
 
     for (int i = 0; i < fringe->x->n; i++) {
@@ -150 +165 @@
         psImage *subImage = psImageSubset(image, region);
         psImage *subMask = psImageSubset(mask, region);
-        psImageStats(stats, subImage, subMask, maskVal);
-
-        fPt[i] = stats->robustMedian;
-        dfPt[i] = stats->robustStdev;
-
-        psTrace(__func__, 7, "[%d:%d,%d:%d]: %f %f\n", region.x0, region.x1, region.y0, region.y1,
-                fPt[i], dfPt[i]);
+        psImageStats(medianSd, subImage, subMask, maskVal);
+        psFree(subImage);
+        psFree(subMask);
+
+        int xSky = xPt[i] / (float)image->numCols * (float)sky->numCols;
+        int ySky = yPt[i] / (float)image->numRows * (float)sky->numRows;
+
+        fPt[i] = medianSd->sampleMedian - sky->data.F32[ySky][xSky];
+        dfPt[i] = 1.0 / medianSd->sampleStdev;
+
+        psTrace(__func__, 7, "[%d:%d,%d:%d]: %f %f\n", (int)region.x0, (int)region.x1,
+                (int)region.y0, (int)region.y1, fPt[i], dfPt[i]);
     }
 
@@ -162 +182 @@
 }
 
-// XXX include the fringe error (fringe->df) in the fit?
-pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes)
-{
-    pmFringeScale *scale = pmFringeScaleAlloc(fringes->n); // The resultant fringe scales
-
-    int numCoeffs = fringes->n + 1;
-    int numPoints = science->f->n;
-
-    psImage *A = psImageAlloc(numCoeffs, numCoeffs, PS_TYPE_F64);
-    psVector *B = psVectorAlloc(numCoeffs, PS_TYPE_F64);
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// pmFringeScale
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+static void fringeScaleFree(pmFringeScale *scale)
+{
+    psFree(scale->coeff);
+    psFree(scale->coeffErr);
+    return;
+}
+
+pmFringeScale *pmFringeScaleAlloc(int nFringeFrames)
+{
+    pmFringeScale *scale = psAlloc(sizeof(pmFringeScale));
+    (void)psMemSetDeallocator(scale, (psFreeFunc)fringeScaleFree);
+
+    scale->nFringeFrames = nFringeFrames;
+    scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
+    scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
+    scale->coeff->n = nFringeFrames + 1;
+    scale->coeff->n = nFringeFrames + 1;
+
+    return scale;
+}
+
+// Determine the fringe scales through solving the least-squares problem
+static bool scaleMeasure(pmFringeScale *scale, // Scale to return
+                         pmFringeStats *science, // The fringe measurements for the science image
+                         psArray *fringes // Array of fringe measurements for the templates
+                        )
+{
+    assert(scale);
+    assert(science);
+    assert(fringes);
+    assert(scale->nFringeFrames == fringes->n);
+
+    psVector *mask = science->regions->mask; // The region mask
+
+    int numCoeffs = fringes->n + 1;     // Number of coefficients: scales for the templates plus a background
+    int numPoints = science->regions->nRequested; // Number of points (i.e., fringe measurements)
+
+    psImage *A = psImageAlloc(numCoeffs, numCoeffs, PS_TYPE_F64); // The least-squares matrix
+    psVector *B = psVectorAlloc(numCoeffs, PS_TYPE_F64); // The least-squares vector
     B->n = numCoeffs;
 
+    // Generate the least-squares matrix and vector
     for (int i = 0; i < numCoeffs; i++) {
-        psVector *fringe1 = NULL;
+        psVector *fringe1 = NULL;       // A fringe measurement
         if (i != 0) {
             pmFringeStats *fringe = fringes->data[i - 1];
             fringe1 = fringe->f;
         }
-        for (int j = 0; j < numCoeffs; j++) {
-            psVector *fringe2 = NULL;
+
+        // Fill in the upper part of the matrix
+        for (int j = i; j < numCoeffs; j++) {
+            psVector *fringe2 = NULL;   // Another fringe measurement
             if (j != 0) {
                 pmFringeStats *fringe = fringes->data[j - 1];
@@ -187 +243 @@
             }
 
-            double sum = 0.0;
+            double matrix = 0.0;        // The matrix sum
             for (int k = 0; k < numPoints; k++) {
-                psF32 f1 = (fringe1 == NULL) ? 1.0 : fringe1->data.F32[k];
-                psF32 f2 = (fringe2 == NULL) ? 1.0 : fringe2->data.F32[k];
-                sum += f1*f2;
-            }
-            A->data.F64[i][j] = sum;
-        }
-
-        double sum = 0.0;
-        for (int j = 0; j < numPoints; j++) {
-            psF32 f1 = (fringe1 == NULL) ? 1.0 : fringe1->data.F32[j];
-            psF32 f2 = science->f->data.F32[j];
-            sum += f1*f2;
-        }
-        B->data.F64[i] = sum;
-    }
-
+                if (!mask->data.U8[k]) {
+                    psF32 f1 = (fringe1) ? fringe1->data.F32[k] : 1.0; // Contribution from i fringe
+                    psF32 f2 = (fringe2) ? fringe2->data.F32[k] : 1.0; // Contribution from j fringe
+                    psF32 dsInv = science->df->data.F32[k]; // 1 / sigma
+                    matrix += f1 * f2 * dsInv * dsInv;
+                }
+            }
+            A->data.F64[i][j] = matrix;
+        }
+
+        // Use symmetry to fill in the lower part of the matrix
+        for (int j = 0; j < i; j++) {
+            A->data.F64[i][j] = A->data.F64[j][i];
+        }
+
+        double vector = 0.0;            // The vector sum
+        for (int k = 0; k < numPoints; k++) {
+            if (!mask->data.U8[k]) {
+                psF32 f1 = (fringe1) ? fringe1->data.F32[k] : 1.0; // Contribution from fringe 1
+                psF32 s = science->f->data.F32[k]; // Contribution from science measurement
+                psF32 dsInv = science->df->data.F32[k]; // 1 / sigma
+                vector += f1 * s * dsInv * dsInv;
+            }
+        }
+        B->data.F64[i] = vector;
+    }
+
+    if (psTraceGetLevel(__func__) >= 5) {
+        printf("From %d points:\n", numPoints);
+        for (int i = 0; i < numCoeffs; i++) {
+            for (int j = 0; j < numCoeffs; j++) {
+                printf("%.2e ", A->data.F64[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    // Solve the least-squares equation
     if (! psGaussJordan(A, B)) {
         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
-    }
-
+        return false;
+    }
+
+    // Copy the results over
     for (int i = 0; i < numCoeffs; i++) {
         scale->coeff->data.F32[i] = B->data.F64[i];
@@ -214 +294 @@
     }
 
+    return true;
+}
+
+// Measure the fringe differences for each region
+static bool fringeScaleDiffs(psVector *diff, // Vector of differences
+                             pmFringeStats *science, // Science fringe measurements
+                             psArray *fringes, // Template fringe measurements
+                             pmFringeScale *scale // Fringe scales
+                            )
+{
+    assert(diff);
+    assert(diff->type.type == PS_TYPE_F32);
+    assert(science);
+    assert(fringes);
+    assert(scale);
+    assert(diff->n == science->regions->nRequested);
+    assert(fringes->n == scale->nFringeFrames);
+
+    psVector *mask = science->regions->mask; // The region mask
+
+    for (int i = 0; i < diff->n; i++) {
+        if (!mask->data.U8[i]) {
+            float difference = science->f->data.F32[i] - scale->coeff->data.F32[0];
+            for (int j = 0; j < fringes->n; j++) {
+                pmFringeStats *fringe = fringes->data[j]; // The fringe of interest
+                difference -= scale->coeff->data.F32[j + 1] * fringe->f->data.F32[i];
+            }
+            diff->data.F32[i] = difference * difference * science->df->data.F32[i] * science->df->data.F32[i];
+        }
+    }
+
+    return true;
+}
+
+// Clip regions based on the differences; return the number masked
+static int clipRegions(psVector *diffs, // Differences
+                       psVector *mask,  // Region mask
+                       float rej        // Rejection limit in standard deviations
+                      )
+{
+    assert(diffs);
+    assert(diffs->type.type == PS_TYPE_F32);
+    assert(mask);
+    assert(mask->type.type == PS_TYPE_U8);
+    assert(diffs->n == mask->n);
+
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics
+    psVectorStats(stats, diffs, NULL, mask, 1);
+    float middle = stats->sampleMedian; // The middle of the distribution
+    float thresh = rej * 0.74 * (stats->sampleUQ - stats->sampleLQ); // The rejection threshold
+    psFree(stats);
+
+    int numClipped = 0;                 // Number clipped
+    for (int i = 0; i < diffs->n; i++) {
+        psTrace(__func__, 10, "Region %d (%d): %f\n", i, mask->data.U8[i], diffs->data.F32[i]);
+        if (!mask->data.U8[i] && fabs(diffs->data.F32[i]) > middle + thresh) {
+            psTrace(__func__, 5, "Masking %d: %f\n", i, diffs->data.F32[i]);
+            mask->data.U8[i] = 1;
+            numClipped++;
+        }
+    }
+
+    return numClipped;
+}
+
+
+// XXX include the fringe error (fringe->df) in the fit?
+pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes, float rej,
+                                    unsigned int nIter, float keepFrac)
+{
+    pmFringeRegions *regions = science->regions; // The fringe regions
+    int numRegions = regions->nRequested; // Number of regions
+    if (!regions->mask) {
+        regions->mask = psVectorAlloc(numRegions, PS_TYPE_U8);
+        regions->mask->n = numRegions;
+        psVectorInit(regions->mask, 0);
+    }
+
+    psVector *mask = regions->mask;     // The region mask
+    psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Median statistics
+    unsigned int numClipped = 0;        // Total number clipped
+    psVector *diff = psVectorAlloc(numRegions, PS_TYPE_F32); // The differences between obs. and pred.
+    diff->n = numRegions;
+
+    pmFringeScale *scale = pmFringeScaleAlloc(fringes->n); // The fringe scales
+
+    // Get rid of bad data points
+    for (int i = 0; i < fringes->n; i++) {
+        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
+        for (int j = 0; j < numRegions; j++) {
+            if (!finite(fringe->f->data.F32[j])) {
+                mask->data.U8[j] = 1;
+                psTrace(__func__, 9, "Masking region %d because not finite in fringe %d.\n", j, i);
+            }
+        }
+    }
+
+    // Get rid of the extreme outliers by assuming most of the points are somewhat clustered
+    psVectorStats(median, science->f, NULL, NULL, 0);
+    scale->coeff->data.F32[0] = median->sampleMedian;
+    for (int i = 0; i < fringes->n; i++) {
+        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
+        psVectorStats(median, fringe->f, NULL, NULL, 0);
+        scale->coeff->data.F32[0] -= median->sampleMedian;
+        scale->coeff->data.F32[i] = 0.0;
+    }
+    psFree(median);
+    fringeScaleDiffs(diff, science, fringes, scale);
+    numClipped = clipRegions(diff, mask, 3.0*rej);
+    psTrace(__func__, 4, "%d regions clipped in initial pass.\n", numClipped);
+
+    unsigned int iter = 0;              // Iteration number
+    unsigned int iterClip = 0;          // Number clipped in this iteration
+    do {
+        iter++;
+        scaleMeasure(scale, science, fringes); // The scales
+        psTrace(__func__, 1, "Fringe scales after iteration %d:\n", iter);
+        psTrace(__func__, 1, "Background: %f %f\n", scale->coeff->data.F32[0], scale->coeffErr->data.F32[0]);
+        for (int i = 0; i < scale->nFringeFrames; i++) {
+            psTrace(__func__, 1, "%d: %f %f\n", i, scale->coeff->data.F32[i + 1],
+                    scale->coeffErr->data.F32[i + 1]);
+        }
+
+        fringeScaleDiffs(diff, science, fringes, scale);
+        iterClip = clipRegions(diff, mask, rej); // Number clipped
+        numClipped += iterClip;
+        psTrace(__func__, 9, "Clipped: %d\tFrac: %f\n", iterClip, (float)numClipped/(float)numRegions);
+    } while (iterClip > 0 && iter < nIter && (float)numClipped/(float)numRegions <= 1.0 - keepFrac);
+    psFree(diff);
+
+    // A final iteration with the last clipping
+    scaleMeasure(scale, science, fringes);
+
     return scale;
 }
 
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Fringe correction
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
 // XXX note that this modifies the input fringe images
-psImage *pmFringeCorrect(pmReadout *readout, pmFringeRegions *fringes, psArray *fringeImages, psArray *fringeStats, psMaskType maskVal)
+psImage *pmFringeCorrect(pmReadout *readout, pmFringeRegions *fringes, psArray *fringeImages,
+                         psArray *fringeStats, psMaskType maskVal, float rej,
+                         unsigned int nIter, float keepFrac)
 {
     // measure the fringe stats for the science frame and solve for the scales
     pmFringeStats *scienceStats = pmFringeStatsMeasure(fringes, readout, maskVal);
-    pmFringeScale *scale = pmFringeScaleMeasure(scienceStats, fringeStats);
+
+    if (psTraceGetLevel(__func__) > 9) {
+        for (int i = 0; i < fringes->nRequested; i++) {
+            printf("%f", scienceStats->f->data.F32[i]);
+            for (int j = 0; j < fringeStats->n; j++) {
+                pmFringeStats *fringe = fringeStats->data[j];
+                printf("\t%f", fringe->f->data.F32[i]);
+            }
+            printf("\n");
+        }
+    }
+
+    pmFringeScale *scale = pmFringeScaleMeasure(scienceStats, fringeStats, rej, nIter, keepFrac);
     psFree(scienceStats);
 
     psTrace(__func__, 7, "Fringe solution:\n");
-    for (int i = 0; i < fringeImages->n; i++) {
+    for (int i = 0; i < fringeImages->n + 1; i++) {
         psTrace(__func__, 7, "%d: %f %f\n", i, scale->coeff->data.F32[i],
                 scale->coeffErr->data.F32[i]);
@@ -234 +465 @@
     // XXX we could save data space by making the first image the output image
     psImage *sumFringe = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
+    //psBinaryOp(sumFringe, sumFringe, "+", psScalarAlloc(scale->coeff->data.F32[0], PS_TYPE_F32));
     for (int i = 0; i < fringeImages->n; i++) {