Changeset 34848

Timestamp:

Dec 18, 2012, 2:12:42 PM (13 years ago)

Author:

watersc1

Message:

Added linear scaling code, which is off by default. These models do look better than with no scaling, but they're not as good as my manual models. I'm not sure what is different, but I wanted to get this into SVN.

File:

• trunk/ppStack/src/ppStackCombineAlternate.c (modified) (3 diffs)

trunk/ppStack/src/ppStackCombineAlternate.c

@@ -23 +23 @@
     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 +87 @@
     pmFPAfileClose(file,view);
   }
+
+
+  // Do combination
   if (!pmStackSimpleMedianCombine(bkgRO,inputs)) {
     psFree(inputs);
@@ -113 +124 @@
   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);
+}
+
+