Changeset 41892 for trunk/psModules/src/detrend/pmPattern.c

Timestamp:

Nov 4, 2021, 6:05:18 PM (5 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-dev-20210817 (pmPattern updates, new inverse transform extra orders api, forward transform uses ORD, pmSourceIO_CMF.c.in conversion to pmFitsTableNew)

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/detrend/pmPattern.c (modified) (10 diffs)

trunk/psModules/src/detrend/pmPattern.c

@@ -9 +9 @@
 #define PATTERN_ROW_BKG_FIX 1
 
+/* some in-line notes:
+
+   patternMaskRow sets the data value to NAN, inconsistent with new plan
+   
+   here is the outline of pmPatternRow
+
+   * at this point, we have already done overscan subtraction, right?
+
+   * measure stats on the full cell (MEDIAN, STDEV)
+   ** subsample?
+   ** if it fails, it masks the entire cell and set the value to NAN
+
+   * calculate an upper and lower threshold (median +/- T * sigma)
+   * define a normalized x-coordinate ('index') : 
+   ** see note below on chebys
+
+   * each row is treated independently
+   * pixels are masked for the fit if they are out-of-range 
+     or if they are already masked
+
+   ** note that the clipping threshold will be larger if there 
+      are pixels which have astronomical structures
+      
+      a possible better option would be to set the threshold based on the median
+      and a sigma calculated from Poisson stats (do we know the gain?)
+      
+   ** fit is allowed to proceed if even N+1 pixels exist, which is clearly too low
+   
+   ** Remaining pixels are fitted with clip-fit 
+
+   ** solution is subtracted from the data
+   (this is implemented with psPolynomial1DEvalVector)
+   perhaps faster if we fixed the order to 2 and hardwired the result
+   
+   * after each row is fitted, the intercept (A value) is fitted
+   as a function of the y-coordinate and the result is subtracted
+
+   * the slope value is also fitted as a function of the column 
+     and added back in -- I'm not sure I understand this step.
+
+   *****************
+
+   ** what we calculate are related to chebychevs (domain is -1 : +1)
+   *** T0(x) = 1
+   *** T1(x) = x
+   *** T2(x) = 2x^2 - 1
+
+   *** we calculate y = A + Bx + Cx^2
+
+   a_0 + a_1 x + a_2 (2x^2 - 1) = A + B + Cx^2
+
+   a_1       = B
+   a_0 - a_2 = A
+   2 a_2     = C
+
+   a_0       = A + C/2
+   a_1       = B
+   a_2       = C/2
+
+   *****************
+   
+   I have 3 goals in re-working the code:
+   
+   1) improve overall speed
+   2) improve reliability of the fit
+   3) skip fit if we can
+
+   Let's assume the signal in the cell is light + bias drift
+
+   The bias drift has an amplitude of ~5 - 10 DN
+
+   That makes a detectable source with ~N * a few counts (multiple pixels in a row)
+   
+   So, the effective flux is ~10 * 5 = 50 DN
+   for which sky level is this value - 3 sigma?
+   
+   50 / sqrt(sky sigma^2 * effective area)
+
+   area ~ 5pixels, sky sigma^2 = sky
+
+   10 * Npix / sqrt(sky * Npix) = 3
+
+   Npix = sky * (S/N)^2 / (peak^2)
+   
+   sky = Npix * peak^2 / SN^2
+
+   if (sky < Npix * peak^2 / SN^2), we should skip:
+   
+   Npix ~ 5
+   peak ~ 10
+   SN ~ 3 (or even less)
+
+   sky < 5 * 100 / 9 = 55 or so
+
+
+   To address these in order:
+   
+   1) speed: 
+   * the analysis threaded is not threaded: thread across cells
+   
+   2) 
+
+ */
 
 // Mask a row as bad
@@ -33 +136 @@
 }
 
+// Comparison and swap functions for sorting values directly
+#define SORT_COMPARE(A,B) (sampleArray[A] < sampleArray[B])
+#define SORT_SWAP(TYPE,A,B) {                   \
+    if (A != B) { \
+        TYPE temp = sampleArray[A];                     \
+        sampleArray[A] = sampleArray[B]; \
+        sampleArray[B] = temp; \
+    } \
+}
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Measurement and application
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                          psStatsOptions clipMean, psStatsOptions clipStdev,
+                          psImageMaskType maskVal, psImageMaskType maskBad);
+
+
+bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                        psStatsOptions clipMean, psStatsOptions clipStdev,
+                        psImageMaskType maskVal, psImageMaskType maskBad);
+
+
+// XXX allow user choice of binned vs unbinned analysis?
 bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh,
+                  psStatsOptions clipMean, psStatsOptions clipStdev,
+                  psImageMaskType maskVal, psImageMaskType maskBad) {
+
+  bool status = false;
+  if (true) {
+    status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
+  } else {
+    status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
+  }
+  return status;
+}
+
+// USE_BACKGROUND_STDEV: if TRUE, the analysis will use the measured robust stdev to clip the out-of-range pixles
+// if FALSE, the stdev will be estimated based on the Poisson statistics of the median (sky level).
+# define USE_BACKGROUND_STDEV 0
+
+bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh,
                   psStatsOptions clipMean, psStatsOptions clipStdev,
                   psImageMaskType maskVal, psImageMaskType maskBad)
@@ -48 +189 @@
     PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
 
+    bool mdok;                          // Status of MD lookup
+
+    pmCell *cell = ro->parent;
+    pmChip *chip = cell->parent;
+    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
+    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
+
     psImage *image = ro->image;         // Image to correct
     psImage *mask = ro->mask;           // Mask for image
@@ -55 +203 @@
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) {
-        psWarning("Unable to calculate statistics on readout; masking entire readout.");
+        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
         psErrorClear();
         psFree(stats);
         psFree(rng);
-        psImageInit(image, NAN);
-        if (mask) {
-            psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK));
-        }
-        if (ro->variance) {
-            psImageInit(image, NAN);
-        }
         return true;
     }
+
+# if (USE_BACKGROUND_STDEV) 
     float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
     float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data
     float background = stats->robustMedian;
+# else
+    // the signal we are looking for is a small variation on top of the background.  if
+    // the background is uniform with only read noise + sky noise, then the pixel-to-pixel
+    // stdev should only be due to known noise sources and predictable.  If the
+    // pixel-to-pixel variations are from other features, then those variations will
+    // probably dominate the row-by-row bias variations.
+
+    // instead of using the image pixel statistics to measure the stdev, lets assume only
+    // dark noise plus poisson sky noise.  we are not carrying in the read noise, but it is
+    // fairly modest for GPC1 (~10 DN)
+
+    // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have
+    // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17
+
+    // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2)
+    // so smaller by a factor of 1.4 than what we predict, which is not very large
+
+    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
+    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
+    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    // If we cannot determine the nominal amplitude, we fall-back on the worst case
+
+    // XXX retrieve noise and gain from 'concepts' if possible
+# define READNOISE 10 /* arbitrary number */
+    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); 
+    float lower = stats->robustMedian - delta; // Lower bound for data
+    float upper = stats->robustMedian + delta; // Upper bound for data
+    float background = stats->robustMedian;
+
+    // if the noise from the background is too large 
+    float significance = nominalAmplitude / sigma;
+   
+    // XXX EAM : arbitrary number
+    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      return true; 
+    }
+    fprintf (stderr, "correcting pattern row background %s: %f - %f - %f : %f %f %f\n", cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+# endif    
+    
     psFree(stats);
     psFree(rng);
@@ -109 +294 @@
     psVectorInit(yaxisMask, 0);
 #endif
+
+    // we really need more than order + 1 points (= 4).
+    // this should be tunable, but let's try 5 - 10%
+    int validNmin = numCols * 0.1;
+
     for (int y = 0; y < numRows; y++) {
         psVectorInit(clipMask, 0);
         data = psImageRow(data, image, y);
         int num = 0;                    // Number of good pixels
+
+        // if the unmasked pixels only span a small range in x then we cannot fit the
+        // 2nd order polynomial variations very well.  Require a minimum fractional range
+        float validXmin = +1;
+        float validXmax = -1;
+
+        // XXX can we do just as well fitting 1/3 of the pixels? (NOT REALLY)
+        // (x % 3) ||
         for (int x = 0; x < numCols; x++) {
-            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
-                data->data.F32[x] < lower || data->data.F32[x] > upper) {
-                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF;
+            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                data->data.F32[x] < lower || data->data.F32[x] > upper) {
+                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF;
             } else {
                 clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0;
                 num++;
+                validXmin = PS_MIN(indices->data.F32[x], validXmin);
+                validXmax = PS_MAX(indices->data.F32[x], validXmax);
             }
         }
-        if (num < order + 1) {
+
+        // XXX how much time is spent in the fitting
+        if (num < validNmin) {
             // Not enough points to fit
             patternMaskRow(ro, y, maskBad);
@@ -131 +333 @@
             continue;
         }
+        // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise?
         if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, data, NULL, indices)) {
             psWarning("Unable to fit polynomial to row %d", y);
@@ -215 +418 @@
           for (int x = 0; x < numCols; x++) {
             image->data.F32[y][x] += solution->data.F32[y];
-            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]);
+            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]);
           }
           corr->data.F64[y][0]  -= solution->data.F32[y];
@@ -241 +444 @@
           for (int x = 0; x < numCols; x++) {
             image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x];
-            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],indices->data.F32[x]);
+            // XXX EAM : this was [x] which is wrong
+            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]);
           }
           corr->data.F64[y][1]  -= solution->data.F32[y] ;
@@ -262 +466 @@
 
     psFree(indices);
+    psFree(clip);
+    psFree(clipMask);
+    psFree(poly);
+    psFree(data);
+
+    return true;
+}
+
+# define NPIX 15
+
+// bin by NPIX in the x-direction to reduce the number of calculations needed to measure
+// the pattern
+bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                  psStatsOptions clipMean, psStatsOptions clipStdev,
+                  psImageMaskType maskVal, psImageMaskType maskBad)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PM_ASSERT_READOUT_IMAGE(ro, false);
+    PS_ASSERT_INT_NONNEGATIVE(order, false);
+    PS_ASSERT_INT_NONNEGATIVE(iter, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
+
+    bool mdok;                          // Status of MD lookup
+
+    pmCell *cell = ro->parent;
+    pmChip *chip = cell->parent;
+    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
+    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
+
+    psImage *image = ro->image;         // Image to correct
+    psImage *mask = ro->mask;           // Mask for image
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+    if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) {
+        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
+        psErrorClear();
+        psFree(stats);
+        psFree(rng);
+        
+        // EAM 20211011 : we used to mask cells which fail the above, but this seems excessive
+        // psImageInit(image, NAN);
+        // if (mask) {
+        //     psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK));
+        // }
+        // if (ro->variance) {
+        //     psImageInit(image, NAN);
+        // }
+        return true;
+    }
+
+    // if USE_BACKGROUND_STDEV is TRUE, the observed standard deviation is used to set the
+    // thresholds.  this is going to be an overestimate if there is any structure in the
+    // image.  If FALSE, the thresholds are set based on poisson stats for the background
+    // level.  We assume the gain is 1, so this is an overestimate if the gain is > 1
+
+# if (USE_BACKGROUND_STDEV) 
+    float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
+    float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data
+    float background = stats->robustMedian;
+# else
+    // the signal we are looking for is a small variation on top of the background.  if
+    // the background is uniform with only read noise + sky noise, then the pixel-to-pixel
+    // stdev should only be due to known noise sources and predictable.  If the
+    // pixel-to-pixel variations are from other features, then those variations will
+    // probably dominate the row-by-row bias variations.
+
+    // instead of using the image pixel statistics to measure the stdev, lets assume only
+    // dark noise plus poisson sky noise.  we are not carrying in the read noise, but it is
+    // fairly modest for GPC1 (~10 DN)
+
+    // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have
+    // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17
+
+    // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2)
+    // so smaller by a factor of 1.4 than what we predict, which is not very large
+
+    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
+    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
+    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    if (!isfinite(nominalAmplitude)) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    // If we cannot determine the nominal amplitude, we fall-back on the worst case
+
+    // retrieve noise and gain from 'concepts' if possible
+# define READNOISE 10 /* arbitrary number */
+    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); 
+    float lower = stats->robustMedian - delta; // Lower bound for data
+    float upper = stats->robustMedian + delta; // Upper bound for data
+    float background = stats->robustMedian;
+
+    // if the noise from the background is too large 
+    float significance = nominalAmplitude / sigma;
+   
+    // XXX EAM : arbitrary number
+    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      psFree(stats);
+      psFree(rng);
+      return true; 
+    }
+    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+# endif    
+
+    psFree(stats);
+    psFree(rng);
+
+    // the vector 'indices' maps the x-coordinate to a range [-1:1].  the element number (i) of indices
+    // related to the x-coordinate (column number) by x = (i + 0.5) * NPIX
+
+    int nSamples = numCols / NPIX;
+
+    psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for fit solutions
+    psVector *xFit = psVectorAlloc(nSamples, PS_TYPE_F32); // x-coordinate for fitting
+    psVector *yFit = psVectorAlloc(nSamples, PS_TYPE_F32); // flux values for fitting
+
+    // xFit elements run from 0 - nSamples, element 'sample' corresponds to the middle of the bin sample*NPIX + 0.5*NPIX
+
+    float norm = 2.0 / (float)numCols;  // Normalisation for indices
+    for (int sample = 0; sample < nSamples; sample ++) {
+        int x = (sample + 0.5)*NPIX;
+        xFit->data.F32[sample] = x * norm - 1.0;
+    }
+    for (int x = 0; x < numCols; x ++) {
+        indices->data.F32[x] = x * norm - 1.0;
+    }
+
+    psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics
+    clip->clipIter = iter;
+    clip->clipSigma = rej;
+    psVector *clipMask = psVectorAlloc(nSamples, PS_TYPE_VECTOR_MASK); // Mask for clipping
+    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit
+    psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit
+
+    psImage *corr = psImageAlloc(order + 1, numRows, PS_TYPE_F64); // Corrections applied
+    psImageInit(corr, NAN);
+
+    // CZW: 2011-11-30
+    // Define the vectors to hold the "x" and "y" slope trends.
+    // Briefly, the slope trend in the y-axis is a due to variations in the 0-th order term
+    // of the PATTERN.ROW fit between individual rows across the cell.  Similarly, the 1-st
+    // order term of the PATTERN.ROW fit defines the trend in the x-axis (as that's what we
+    // are fitting with PATTERN.ROW in the first place).  However, the thing we're trying to
+    // fix with PATTERN.ROW is the detector level bias wiggles.  These should be overlaid on
+    // the true sky level.  Therefore, simply applying the PATTERN.ROW correction will
+    // introduce cell-to-cell sky variations as these two trends are removed.  To avoid this,
+    // We store the 0th and 1st order values used for each row, and then fit a polynomial to
+    // these results.  By re-adding these systematic trends back, we can remove the row-to-row
+    // variations without improperly removing the real sky trend.
+    psVector *yaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the constant term
+    psVector *yaxisMask = psVectorAlloc(numRows, PS_TYPE_VECTOR_MASK); // Mask for rows with no fit
+    psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term
+    psVectorInit(yaxisMask, 0);
+
+    // validNmin is the minimum number of samples needed to measure the trend.
+    // this should be tunable, but let's try 5 - 10%
+    int validNmin = PS_MAX (nSamples * 0.1, order + 2);
+
+    for (int y = 0; y < numRows; y++) {
+        psVectorInit(clipMask, 0);
+        data = psImageRow(data, image, y);
+        int num = 0;                    // Number of good pixels
+
+        // if the unmasked pixels only span a small range in x then we cannot fit the
+        // 2nd order polynomial variations very well.  Require a minimum fractional range
+        float validXmin = +1;
+        float validXmax = -1;
+
+        for (int sample = 0; sample < nSamples; sample ++) {
+
+            // store valid samples in the array to be sorted
+            float sampleArray[NPIX];
+            int seq = 0;
+            for (int j = 0; j < NPIX; j++) {
+                int pix = sample  * NPIX + j; // real pixel elements in x-dir
+                psAssert (pix >= 0, "invalid pix value");
+                psAssert (pix < numCols, "invalid pix value");
+                if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][pix] & maskVal)) continue;
+                if (data->data.F32[pix] < lower || data->data.F32[pix] > upper) continue;
+                sampleArray[seq] = data->data.F32[pix]; // store the value to be sorted
+                seq ++;
+            } 
+            if (seq < 1) {
+                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0xFF;
+                yFit->data.F32[sample] = NAN;
+                continue;
+            }
+            // note that we are treating the x-coordinate as the center
+            // of the binned pixel group, even if some or most pixels have
+            // been masked.  compared to the amplitude of the slope, this
+            // error is small
+            clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0;
+            validXmin = PS_MIN(xFit->data.F32[sample], validXmin);
+            validXmax = PS_MAX(xFit->data.F32[sample], validXmax);
+            num++;
+
+            // PSSORT operates on sampleArray (see define of macro SORT_SWAP above)
+            PSSORT (seq, SORT_COMPARE, SORT_SWAP, float);
+
+            int midPt = 0.5 * seq;
+            if (seq % 2 == 1) { psAssert (midPt >= 0, "invalid midPt"); }
+            if (seq % 2 == 0) { psAssert (midPt >= 1, "invalid midPt"); }
+            psAssert (midPt < NPIX, "invalid midPt");
+
+            float medValue = (seq % 2) ? sampleArray[midPt] : 0.5*(sampleArray[midPt] + sampleArray[midPt-1]);
+            yFit->data.F32[sample] = medValue;
+        }
+
+        // If not enough points are valid, skip
+        if (num < validNmin) {
+            // Ignore this row in our subsequent fits, because the fit failed.
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+            continue;
+        }
+        // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise?
+        if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, yFit, NULL, xFit)) {
+            psWarning("Unable to fit polynomial to row %d", y);
+            psErrorClear();
+            // Ignore this row in our subsequent fits, because the fit failed.
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+            continue;
+        }
+        // Store the results we found for this row.
+        yaxisData->data.F32[y] = poly->coeff[0];
+        xaxisData->data.F32[y] = poly->coeff[1];
+        psTrace("pattern",1,"%d %g %g\n",y,poly->coeff[0],poly->coeff[1]);
+
+        memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector
+        if (!solution) {
+            psWarning("Unable to evaluate polynomial for row %d", y);
+            psErrorClear();
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+            continue;
+        }
+
+        psAssert (solution->n == numCols, "oops");
+        for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] -= solution->data.F32[x];
+            psTrace("pattern",5,"A: %d %d %g\n",x,y,solution->data.F32[x]);
+        }
+        psFree(solution);
+    }
+
+    // Put the global trends back that were removed by the PATTERN.ROW correction.
+    // Set up the indices for the polynomial
+    psVector *yaxisIndices = psVectorAlloc(numRows, PS_TYPE_F32);
+    norm = 2.0 / (float)numRows;
+    for (int y = 0; y < numRows; y++) {
+      yaxisIndices->data.F32[y] = y * norm - 1.0;
+      psTrace("psModules.detrend.pattern",10,"%d %f %f\n",y,yaxisIndices->data.F32[y],yaxisData->data.F32[y]);
+    }
+
+    // Fit the trend of the constant term, producing the y-axis global trend
+    psStatsInit(clip);
+    psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
+    if (!psVectorClipFitPolynomial1D(yaxisPoly,clip,yaxisMask,0xFF,yaxisData, NULL, yaxisIndices)) {
+      psWarning("Unable to fit polynomial to y-axis trend");
+      psErrorClear();
+      // If we've failed, we need to do something, so add back in the background level, and
+      // expect that the final image will have background mismatches.
+      for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+          image->data.F32[y][x] += background;
+        }
+        corr->data.F64[y][0]  -= background;
+      }
+    } else {
+      psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices);
+      if (!solution) {
+        psWarning("Unable to evaluate polynomial");
+        psErrorClear();
+        // If we've failed, we need to do something, so add back in the background level, and
+        // expect that the final image will have background mismatches.
+        for (int y = 0; y < numRows; y++) {
+          for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] += background;
+          }
+          corr->data.F64[y][0]  -= background;
+        }
+      } else {
+        psAssert (solution->n == numRows, "oops");
+        for (int y = 0; y < numRows; y++) {
+          for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] += solution->data.F32[y];
+            // XXX EAM : this was [x], which is wrong
+            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]);
+          }
+          corr->data.F64[y][0]  -= solution->data.F32[y];
+        }
+      }
+      psFree(solution);
+    }      
+
+    // Fit the trend of the linear term, producing the x-axis global trend
+    // We can use the same mask vector, as the same rows failed the row-fit earlier.
+    psStatsInit(clip);
+    psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
+    if (!psVectorClipFitPolynomial1D(xaxisPoly,clip,yaxisMask,0xFF,xaxisData, NULL, yaxisIndices)) {
+      psWarning("Unable to fit polynomial to x-axis trend");
+      psErrorClear();
+    } else {
+      psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices);
+      if (!solution) {
+        psWarning("Unable to evaluate polynomial");
+        psErrorClear();
+      } else {
+        psAssert (solution->n == numRows, "oops");
+        for (int y = 0; y < numRows; y++) {
+          for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x];
+            // XXX EAM : this was set to [x] which is wrong (numCols > numRows)
+            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]);
+          }
+          corr->data.F64[y][1]  -= solution->data.F32[y] ;
+        }
+      }
+      psFree(solution);
+    }
+    psFree(yaxisPoly);
+    psFree(xaxisPoly);
+    psFree(yaxisIndices);
+    psFree(yaxisMask);
+    psFree(yaxisData);
+    psFree(xaxisData);
+    
+    psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE,
+                       "Pattern row correction", corr);
+    psFree(corr);
+
+    psFree(indices);
+    psFree(xFit);
+    psFree(yFit);
     psFree(clip);
     psFree(clipMask);
@@ -1316 +1855 @@
     return true;
 }
+