Index: branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c
===================================================================
--- branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c	(revision 41808)
+++ branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c	(revision 41811)
@@ -134,4 +134,14 @@
     }
     return;
+}
+
+// Comparison and swap functions for sorting values directly
+#define SORT_COMPARE(A,B) (value[A] < value[B])
+#define SORT_SWAP(A,B) { \
+    if (A != B) { \
+        float temp = value[A]; \
+        value[A] = value[B]; \
+        value[B] = temp; \
+    } \
 }
 
@@ -313,4 +323,303 @@
 	
 #endif
+        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();
+            patternMaskRow(ro, y, maskBad);
+#ifdef PATTERN_ROW_BKG_FIX
+	    yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+#endif
+            continue;
+        }
+
+        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);
+    }
+
+#ifdef PATTERN_ROW_BKG_FIX
+    // 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 {
+	for (int y = 0; y < numRows; y++) {
+	  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]);
+	  }
+	  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 {
+	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];
+	    psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],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);
+    // End PATTERN_ROW_BKG_FIX global trend replacement
+#endif 
+    
+    psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE,
+                       "Pattern row correction", corr);
+    psFree(corr);
+
+    psFree(indices);
+    psFree(clip);
+    psFree(clipMask);
+    psFree(poly);
+    psFree(data);
+
+    return true;
+}
+
+// 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);
+
+    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; masking entire readout.");
+        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
+
+    # define READNOISE 10
+    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    float lower = stats->robustMedian - thresh * sigma; // Lower bound for data
+    float upper = stats->robustMedian + thresh * sigma; // Upper bound for data
+    float background = stats->robustMedian;
+# endif    
+    
+    pmCell *cell = ro->parent;
+    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
+    fprintf (stderr, "pattern row background %s: %f - %f - %f\n", cellName,lower, background, upper);
+
+
+    // XXX add code to skip fit if background stdev is large compared to expected signal
+    // XXX add a constraint to the fit so the |amplitude| < x (x ~ 20 - 25)
+
+    psFree(stats);
+    psFree(rng);
+
+# define NPIX 15    
+
+    // Indices are distributed [-1:1] [-1 = 0, +1 = numCols = indices[nSamples]
+    int nSamples = numCols / NPIX;
+    psVector *indices = psVectorAlloc(nSamples, PS_TYPE_F32); // Indices for fitting
+    psVector *fitData = psVectorAlloc(nSAmples, PS_TYPE_F32); // Data to fit
+
+    // indices 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;
+        indices->data.F32[sample] = x * norm - 1.0;
+    }
+
+    psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics
+    // XXX clip->clipIter = iter;
+    clip->clipIter = 1; // XXX skip iteration for a test
+    clip->clipSigma = rej;
+    psVector *clipMask = psVectorAlloc(numCols, 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);
+
+#ifdef PATTERN_ROW_BKG_FIX
+    // 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);
+#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 bookkeeping might be easier if this loop is over elements of 'indices'
+	// XXX can we do just as well fitting 1/3 of the pixels? (NOT REALLY)
+	// (x % 3) ||
+        for (int sample = 0; sample < nSamples; sample ++) {
+
+	    // store valid samples in the array
+	    float sampleArray[NPIX];
+	    int seq = 0;
+	    for (int j = 0; j < NPIX; j++) {
+		int pix = sample  * NPIX + j;
+		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];
+		seq ++;
+            } 
+	    if (seq < 1) {
+		clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0xFF;
+            } else {
+                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0;
+                num++;
+		validXmin = PS_MIN(indices->data.F32[sample], validXmin);
+		validXmax = PS_MAX(indices->data.F32[sample], validXmax);
+            }
+
+	    PSSORT (seq, sampleArray
+
+        }
+
+	// XXX how much time is spent in the fitting
+        if (num < validNmin) {
+            // Not enough points to fit
+            patternMaskRow(ro, y, maskBad);
+	    // 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, data, NULL, indices)) {
+            psWarning("Unable to fit polynomial to row %d", y);
+            psErrorClear();
+            patternMaskRow(ro, y, maskBad);
+	    // 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
@@ -1476,2 +1785,3 @@
     return true;
 }
+