Changeset 34800 for trunk/psModules

Timestamp:

Dec 11, 2012, 2:04:31 PM (13 years ago)

Author:

watersc1

Message:

Merge from my branch including background restoration and stack stage projection cell binned images.

Location:

trunk

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/camera/pmFPABin.c (modified) (4 diffs)
• psModules/src/camera/pmFPAfileIO.c (modified) (1 diff)
• psModules/src/detrend/pmDark.c (modified) (3 diffs)
• psModules/src/imcombine/pmStack.c (modified) (1 diff)
• psModules/src/imcombine/pmStack.h (modified) (1 diff)
• psModules/src/objects/Makefile.am (modified) (1 diff)
• psModules/src/objects/pmFootprintCullPeaks.c (modified) (1 diff)

trunk/psModules/src/camera/pmFPABin.c

@@ -43 +43 @@
     }
 
+    int Nbits = (int) (ceil(log(maskVal)/log(2)) + 1);
+    int *bitcounter = malloc(sizeof(int) * Nbits);
+    int pxlcount;
+
     int xLast = numColsIn - 1, yLast = numRowsIn - 1; // Last index
     int yStart = psImageBinningGetFineY(binning, 0); // Starting input y for binning
@@ -55 +59 @@
             float sum = 0.0;            // Sum of pixels
             int numPix = 0;             // Number of pixels
+
+            for (int j = 0; j < Nbits; j++) { // Reset bit counter
+              bitcounter[j] = 0;
+            }
+            pxlcount = 0;
+            
             for (int y = yStart; y < yStop; y++) {
                 for (int x = xStart; x < xStop; x++) {
+                  if (inMask && (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] != 0)) {
+                      for (int j = 0; j < Nbits; j++) {
+                        psImageMaskType M = (psImageMaskType) pow(2,j);
+                        if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & M) {
+                          bitcounter[j]++;
+                        }
+                      }
+                    }
+                  
                     if (inMask && (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
                         continue;
@@ -65 +84 @@
                     sum += inImage->data.F32[y][x];
                     numPix++;
+
+
                 }
             }
-
+            
+            
             // Values to set
             float imageValue;
@@ -79 +101 @@
             }
             outImage->data.F32[yOut][xOut] = imageValue;
-            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
+            //            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
+            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = 0;
+            for (int j = 0; j < Nbits; j++) {
+              if (bitcounter[j] > 0.5 * pxlcount) {
+                outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] |= (int) pow(2,j);
+              }
+            }
             xStart = xStop;
         }

trunk/psModules/src/camera/pmFPAfileIO.c

@@ -76 +76 @@
     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
         pmFPAfile *file = item->data.V;
-
         switch (place) {
           case PM_FPA_BEFORE:

trunk/psModules/src/detrend/pmDark.c

@@ -353 +353 @@
 
     // retrieve the required parameter vectors
-    psArray *values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
-    psAssert(values, "values not supplied");
+    psArray *in_values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
+    psAssert(in_values, "values not supplied");
     psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK");
     psAssert(roMask, "roMask not supplied");
-    psVector *orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
-    psAssert(orders, "orders not supplied");
-
-    psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting
+    psVector *max_orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
+    psAssert(max_orders, "orders not supplied");
+
+    psArray *values_set = psArrayAlloc(max_orders->n);
+    psArray *poly_set = psArrayAlloc(max_orders->n);
+    psVector *logL = psVectorAlloc(max_orders->n,PS_TYPE_F64);
+
+    for (int i = 0; i < max_orders->n; i++) {
+      psVector *orders = psVectorAlloc(i+1,PS_TYPE_U8);
+      for (int j = 0; j < orders->n; j++) {
+        orders->data.U8[j] = max_orders->data.U8[j];
+      }
+      poly_set->data[i] =  psPolynomialMDAlloc(orders); // Polynomial for fitting
+      
+      psArray *values = psArrayAlloc(in_values->n);
+      
+      for (int j = 0; j < values->n; j++) {
+        psVector *these_values = psVectorAlloc(i+1,PS_TYPE_F32);
+        psVector *input_values = in_values->data[j];
+
+        for (int k = 0; k < orders->n; k++) {
+          these_values->data.F32[k] = input_values->data.F32[k];
+        }
+        values->data[j] = these_values;
+      }
+      values_set->data[i] = values;
+    }
 
     // retrieve the norm vector, if supplied
@@ -383 +406 @@
     }
 
-    pmDarkVisualInit(values);
+    pmDarkVisualInit(values_set->data[max_orders->n - 1]);
 
     pmReadout *outReadout = output->readouts->data[0];
@@ -423 +446 @@
             }
 
-            if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
+            int k_best = 0;
+            for (int k = 0; k < max_orders->n; k++) {
+              psPolynomialMD *poly = poly_set->data[k];
+              psArray *values = values_set->data[k];
+              
+              if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
                 psErrorClear();         // Nothing we can do about it
                 psVectorInit(poly->coeff, NAN);
-            }
-
-            pmDarkVisualPixelFit(pixels, mask);
-            pmDarkVisualPixelModel(poly, values);
-
-            for (int k = 0; k < poly->coeff->n; k++) {
+              }
+
+              pmDarkVisualPixelFit(pixels, mask);
+              pmDarkVisualPixelModel(poly, values);
+
+              // Insert math here to choose optimum model.
+              logL->data.F64[k] = 0.0;
+              psPolynomialMD *polySig = poly_set->data[0];
+              for (int m = 0; m < poly->deviations->n; m++) {
+                logL->data.F64[k] += pow(poly->deviations->data.F32[m] / polySig->stdevFit,2);
+/*              if ((xOut == 20) && (yOut == 256)) { */
+/*                psTrace("psModules.detrend",3,"pmDarkCombine DEV: %d %d: input %d models: Norders: %d logL: %g value: %g\n", */
+/*                        xOut,yOut,m,k,logL->data.F64[k],poly->deviations->data.F32[m]); */
+/*              } */
+              }
+              if (k > 0) {
+                if ( ( logL->data.F64[k - 1] - logL->data.F64[k] ) > 1) { // Hard coded criterion for a ~5% limit with one degree of freedom
+                  k_best = k;
+                }
+              }
+              if ((xOut <= 600) && (yOut <= 600)) {
+                psTrace("psModules.detrend",3,"pmDarkCombine: %d %d: models: Norders: %d logL: %g BestOrders: %d\n",
+                        xOut,yOut,k,logL->data.F64[k],k_best);
+              }
+            }
+            if (k_best > 1) {
+              k_best = 1;
+            }
+/*          k_best = 1; */
+            // Select the polynomial that seems best.
+            psPolynomialMD *poly = poly_set->data[k_best];
+              
+            //            for (int k = 0; k < poly->coeff->n; k++) {
+            for (int k = 0; k < max_orders->n + 1; k++) { // There is one more coefficient than is stored here.
                 pmReadout *ro = output->readouts->data[k]; // Readout of interest
-                ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+                if (k < poly->coeff->n) {
+                  ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+                }
+                else {
+                  ro->image->data.F32[yOut][xOut] = 0.0;
+                }
             }
             counts->data.U16[yOut][xOut] = poly->numFit;

trunk/psModules/src/imcombine/pmStack.c

@@ -1358 +1358 @@
 }
 
+bool pmStackSimpleMedianCombine(
+                                pmReadout *combined,
+                                psArray *input) {
+  int num = input->n;
+  //  int numCols, numRows;
+  int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+  int xSize, ySize;                   // Size of the output image
+
+  psArray *stack = psArrayAlloc(num); // Stack of readouts  
+  for (int i = 0; i < num; i++) {
+    //    pmStackData *data = input->data[i]; // Stack data for this input
+    pmReadout *ro = input->data[i]; // data->readout;  // Readout of interest
+    if (!ro) {
+      continue;
+    }
+    stack->data[i] = psMemIncrRefCounter(ro);
+  }    
+
+  if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
+                              stack)) {
+    psError(psErrorCodeLast(), false, "Input stack is not valid.");
+    psFree(stack);
+    return false;
+  }
+
+  psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32);
+  psStats  *stats     = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+
+  for (int y = minInputRows; y < maxInputRows; y++) {
+    for (int x = minInputCols; x < maxInputCols; x++) {
+      for (int i = 0; i < input->n; i++) {
+        pmReadout *ro  = stack->data[i];
+        psImage *image = ro->image;
+        pixelData->data.F32[i] = image->data.F32[y][x];
+      }
+      if (!psVectorStats(stats,pixelData,NULL,NULL,0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate median");
+        psFree(stats);
+        psFree(pixelData);
+        psFree(stack);
+        return(false);
+      }
+      combined->image->data.F32[y][x] = stats->robustMedian;
+      combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+    }
+  }
+  
+  psFree(stats);
+  psFree(pixelData);
+  psFree(stack);
+  return (true);
+}
+
 /// Stack input images
 bool pmStackCombine(

trunk/psModules/src/imcombine/pmStack.h

@@ -41 +41 @@
                               float addVariance ///< Additional variance when rejecting
     );
+/// Stack input images simply
+bool pmStackSimpleMedianCombine(pmReadout *combined, ///< Combined readout (output)
+                                psArray *input       ///< Input array of pmStackData
+                                );
 
 /// Stack input images

trunk/psModules/src/objects/Makefile.am

@@ -111 +111 @@
         pmSourceOutputs.h \
         pmSourceIO.h \
+        pmSourceSatstar.h \ 
         pmSourcePlots.h \
         pmSourceVisual.h \

trunk/psModules/src/objects/pmFootprintCullPeaks.c

@@ -91 +91 @@
 
         // max flux is above threshold for brightest peak
-        pmPeak *maxPeak = fp->peaks->data[0];
+      pmPeak *maxPeak = NULL;
+      for (int i = 0; i < fp->peaks->n; i++) {
+        pmPeak *testPeak = fp->peaks->data[i];
+        float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
+        
+        if (isfinite(this_peak)) {
+          maxPeak = fp->peaks->data[i];
+          break;
+        }
+      }
+      psAssert(maxPeak,"maxPeak was not set in these peaks");
+      //      = fp->peaks->data[0];
         float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux;