Index: trunk/ppStack/src/ppStackCombineAlternate.c
===================================================================
--- trunk/ppStack/src/ppStackCombineAlternate.c	(revision 34846)
+++ trunk/ppStack/src/ppStackCombineAlternate.c	(revision 34848)
@@ -23,8 +23,16 @@
     pmFPAfileClose(file,view);
   }
+
+  // Scaling block
+  if (!ppStackLinearScale(inputs, config)) {
+    psFree(inputs);
+    return(false);
+  }
+
   if (!pmStackSimpleMedianCombine(outRO,inputs)) {
     psFree(inputs);
     return(false);
   }
+
 #if 0
   if (!ppStackWriteImage("/tmp/test_forced.median.fits",
@@ -79,4 +87,7 @@
     pmFPAfileClose(file,view);
   }
+
+
+  // Do combination
   if (!pmStackSimpleMedianCombine(bkgRO,inputs)) {
     psFree(inputs);
@@ -113,3 +124,89 @@
   return(true);
 }
-  
+
+bool ppStackLinearScale (psArray *inputs, pmConfig *config)  {
+  psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE);
+  bool doLinearScaling = psMetadataLookupBool(NULL, recipe, "DO.LINEAR.INPUT.SCALING");
+  if (doLinearScaling) {
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+    
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,inputs)) {
+      psError(psErrorCodeLast(), false, "Input stack is not valid.");
+      psFree(inputs);
+      return false;
+    }
+    
+    // Determine best input.
+    int ref = 0;
+    int Nref = 0;
+    double refScale = 0;
+    
+    for (int i = 0; i < inputs->n; i++) {
+      pmReadout *ro = inputs->data[i];
+      psImage *roImage = ro->image;
+      int Nvalid = 0;
+      double scale = 0;
+
+      for (int y = minInputRows; y < maxInputRows; y++) {
+	for (int x = minInputCols; x < maxInputCols; x++) {
+	  if ((isfinite(roImage->data.F32[y][x]))) {
+	    Nvalid += 1;
+	    scale += roImage->data.F32[y][x];
+	  }
+	}
+      }
+      if ((scale > refScale)&&(Nvalid > 0.9 * Nref)) {
+	ref = i;
+	refScale = scale;
+	Nref = Nvalid;
+      }
+    }
+    fprintf(stderr,"ref: %d %d %f\n",ref,Nref,refScale);
+    // Calculate scaling factors
+    pmReadout *refReadout = inputs->data[ref];
+    psImage *refImage     = refReadout->image;
+    for (int i = 0; i < inputs->n; i++) {
+      pmReadout *ro = inputs->data[i];
+      psImage *roImage = ro->image;
+      double S = 0.0;
+      double Sx = 0.0;
+      double Sy = 0.0;
+      double Sxx = 0.0;
+      double Syy = 0.0;
+      double Sxy = 0.0;
+      double D = 0.0;
+      double offset = 0.0;
+      double scale  = 0.0;
+      
+      for (int y = minInputRows; y < maxInputRows; y++) {
+	for (int x = minInputCols; x < maxInputCols; x++) {
+	  if ((isfinite(refImage->data.F32[y][x]))&&
+	      (isfinite(roImage->data.F32[y][x]))) {
+	    S += 1.0;
+	    Sx += roImage->data.F32[y][x];
+	    Sy += refImage->data.F32[y][x];
+	    
+	    Sxx += pow(roImage->data.F32[y][x],2);
+	    Syy += pow(refImage->data.F32[y][x],2);
+	    Sxy += roImage->data.F32[y][x] * refImage->data.F32[y][x];
+	  }
+	}
+      }
+      
+      D = S * Sxx - Sx * Sx;
+      offset = (Sy * Sxx - Sx * Sxy) / D;
+      scale  = (S * Sxy - Sx * Sy) / D;
+      fprintf(stderr,"Scales: %d %g %g\n",i,offset,scale);
+      // Apply scaling factors
+      for (int y = minInputRows; y < maxInputRows; y++) {
+	for (int x = minInputCols; x < maxInputCols; x++) {
+	  roImage->data.F32[y][x] = offset + scale * roImage->data.F32[y][x];
+	}
+      }
+    }      
+  }
+  return(true);
+}
+
+