Changeset 35149 for branches/czw_branch/20130211/psModules/src/detrend/pmPattern.c

Timestamp:

Feb 11, 2013, 4:50:12 PM (13 years ago)

Author:

watersc1

Message:

Changes that should not conflict with the trunk. This is the camera stage background continuity and the pattern.median example code.

File:

• branches/czw_branch/20130211/psModules/src/detrend/pmPattern.c (modified) (3 diffs)

branches/czw_branch/20130211/psModules/src/detrend/pmPattern.c

@@ -4 +4 @@
 
 #include "pmFPA.h"
-
+#include "pmStack.h"
 #include "pmPattern.h"
 
@@ -1085 +1085 @@
         
       }
-      psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) (%d %d) A: %f B: %f C: %f D: %f",
+      psTrace("psModules.detrend.cont",5, "Continuity OTA: %d (%d %d) (%d %d) A: %f B: %f C: %f D: %f",
               i,x,y,xParity,yParity,
               A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]);         
@@ -1302 +1302 @@
     return true;
 }
+
+
+bool pmPatternMedian(pmChip *chip, const psVector *tweak, psStatsOptions bgStat, psStatsOptions cellStat,
+                     psImageMaskType maskVal, psImageMaskType maskBad, bool doRows, bool doCols)
+{
+  //  fprintf(stderr,"Inside pmPatternMedian!\n");
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_VECTOR_NON_NULL(tweak, false);
+    PS_ASSERT_VECTOR_SIZE(tweak, chip->cells->n, false);
+    PS_ASSERT_VECTOR_TYPE(tweak, PS_TYPE_U8, false);
+    
+    int numCells = tweak->n;            // Number of cells
+    numCells = chip->cells->n;
+    psStats *bgStats      = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics on background
+    psRandom *rng         = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+    psVector *cellMedians = psVectorAlloc(numCells, PS_TYPE_F32); // Mask for means
+    psVector *cellSigmas  = psVectorAlloc(numCells, PS_TYPE_F32);
+    psVector *meanMask    = psVectorAlloc(numCells, PS_TYPE_VECTOR_MASK);
+    // Mask bits
+    enum {
+        PM_PATTERN_IGNORE = 0x01,       // Ignore this cell
+        PM_PATTERN_TWEAK  = 0x02,       // Tweak this cell
+        PM_PATTERN_ERROR  = 0x04,       // Error in calculating background
+        PM_PATTERN_ALL    = 0xFF,       // All causes
+    };
+
+    // Count number of cells to tweak
+    int numTweak = 0;                   // Number of cells to tweak
+    int numIgnore = 0;                  // Number of cells to ignore
+    psStatsInit(bgStats);
+
+    for (int i = 0; i < numCells; i++) { // Calculate cell level statistics, and set ignore mask
+      pmCell *cell = chip->cells->data[i]; // Cell of interest
+
+      if (!cell || !cell->data_exists || !cell->process ||
+          cell->readouts->n == 0 || cell->readouts->n > 1 || !cell->readouts->data[0]) {
+        numIgnore++;
+        meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_IGNORE;
+        cellMedians->data.F32[i] = NAN;
+        cellSigmas->data.F32[i] = NAN;
+        continue;
+      }
+      pmReadout *ro = cell->readouts->data[0];
+      if (tweak->data.U8[i]) {
+        meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_TWEAK;
+        numTweak++;
+      }
+      if (!psImageBackground(bgStats, NULL, ro->image, ro->mask, maskVal, rng)) {
+        psWarning("Unable to measure background for cell %d\n", i);
+        psErrorClear();
+        meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR;
+        continue;
+      }
+      cellMedians->data.F32[i] = bgStats->robustMedian;
+      cellSigmas->data.F32[i]  = bgStats->robustStdev;  
+    }
+    if (numTweak == 0) {
+      //        fprintf(stderr,"Inside pmPatternMedian but I'm claiming nothing to do!\n");
+        // Nothing to do
+        psFree(meanMask);
+        return true;
+    }
+
+    int num = 8;
+    
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE);
+    //    psArray *inputs = psArrayAlloc(num);
+    psVector *inputs    = psVectorAlloc(num,PS_TYPE_S32);
+    psVector *pixelData = psVectorAlloc(inputs->n,PS_TYPE_F32);
+    psVector *pixelMask = psVectorAlloc(inputs->n,PS_TYPE_VECTOR_MASK);
+
+    for (int mode = 0; mode <= 1; mode++) {
+      if ((mode == 0)&&(!doRows)) {
+        continue;
+      }
+      if ((mode == 1)&&(!doCols)) {
+        continue;
+      }
+      for (int i = 0; i < num; i++) { // For each row of cells
+        int k = 0;
+        int numvalid = 0;
+        for (int j = 0; j < numCells; j++) { // Scan the cell list for matches for this row.
+          if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE) {
+            continue;
+          }
+          if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_TWEAK)) {
+            continue;
+          }
+
+          pmCell *cell = chip->cells->data[j];
+          const char *cellName = psMetadataLookupStr(NULL,cell->concepts, "CELL.NAME");
+          int x = cellName[2] - '0';
+          int y = cellName[3] - '0';
+          if (((mode == 0)&&(doRows)&&(y == i))||
+              ((mode == 1)&&(doCols)&&(x == i))) {
+            inputs->data.S32[k] = j; // Add the index to the list of inputs.
+            k++;
+            numvalid++;
+          }
+        }
+        // Allocate the median image and mask
+        pmCell    *exCell  = chip->cells->data[inputs->data.S32[0]];
+        pmReadout *example = exCell->readouts->data[0];
+        pmReadout *median = pmReadoutAlloc(NULL);
+        median->image = psImageAlloc(example->image->numCols,
+                                     example->image->numRows,
+                                     PS_TYPE_F32);
+        median->mask = psImageAlloc(example->image->numCols,
+                                     example->image->numRows,
+                                    PS_TYPE_IMAGE_MASK);
+        // Block to do median.
+        {
+          psImageMaskType badMask = 0x7F;  // Configuable would be better.
+          for (int y = 0 ; y < example->image->numRows; y++) {
+            for (int x = 0; x < example->image->numCols; x++) {
+              int badVal = 0;
+              for (k = 0; k < numvalid; k++) {
+                pmCell *cell   = chip->cells->data[inputs->data.S32[k]];
+                pmReadout *ro  = cell->readouts->data[0];
+                psImage *image = ro->image;
+                psImage *mask  = ro->mask;
+                
+                pixelData->data.F32[k] = image->data.F32[y][x];
+                pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[k] = 0;
+                //#if (0)
+                double cell_median = cellMedians->data.F32[inputs->data.S32[k]];
+                double cell_sigma  = cellSigmas->data.F32[inputs->data.S32[k]];
+                if (fabs(pixelData->data.F32[k] - cell_median) > (3.0 * cell_sigma)) {
+                  pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[k] = 1;
+                  badVal++;
+                }
+                //#endif
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & badMask) {
+                  pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[k] = 1;
+                  badVal++;
+                }
+              }
+              for (k = numvalid; k < num; k++) {
+                pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[k] = 1;
+                badVal++;
+              }
+              if (badVal >= num) {
+                median->image->data.F32[y][x] = NAN;
+              }
+              else {
+                if (!psVectorStats(stats,pixelData,NULL,pixelMask,1)) {
+                  // Die gracefully.
+                  return(false);
+                }
+                int worst_k = -1;
+                double worst_v = 0;
+                for (k = 0; k < num; k++) {
+                  if (pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[k] == 0) {
+                    double v = fabs(stats->sampleMedian - pixelData->data.F32[k]);
+                    if (v > worst_v) {
+                      worst_v = v;
+                      worst_k = k;
+                    }
+                  }
+                }
+                if (worst_k != -1) {
+                  pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[worst_k] = 1;
+                  if (!psVectorStats(stats,pixelData,NULL,pixelMask,1)) {
+                    // Die gracefully.
+                    return(false);
+                  }
+                }
+                median->image->data.F32[y][x] = stats->sampleMedian;
+              }
+            }
+          }
+        } // End median block
+        
+        double correction_median = 0;
+        if (!psImageBackground(bgStats, NULL, median->image, NULL, maskVal, rng)) {
+          psWarning("Unable to measure background for cell section %d\n", i);
+          psErrorClear();
+        }
+        correction_median = bgStats->robustMedian;
+        for (int v = 0; v < median->image->numRows; v++) {
+          for (int u = 0; u < median->image->numCols; u++) {
+            median->image->data.F32[v][u] -= correction_median;
+          }
+        }
+
+        // Block to subtract median
+        for (int j = 0; j < numCells; j++) {
+          if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE) {
+            continue;
+          }
+          if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_TWEAK)) {
+            continue;
+          }
+
+          pmCell *cell = chip->cells->data[j];
+
+          const char *cellName = psMetadataLookupStr(NULL,cell->concepts, "CELL.NAME");
+          int x = cellName[2] - '0';
+          int y = cellName[3] - '0';
+          if (((mode == 0)&&(doRows)&&(y == i))||
+              ((mode == 1)&&(doCols)&&(x == i))) {
+            psLogMsg("psModules.detrend",PS_LOG_DETAIL,"PATTERN.MEDIAN mode: %d Cell Statistics Pre: %d %d %d %s %g %g\n",
+                     mode, j,x,y,cellName,cellMedians->data.F32[j],cellSigmas->data.F32[j]);
+            pmReadout *ro = cell->readouts->data[0];
+            psImage *im = ro->image;
+            for (int v = 0; v < im->numRows; v++) {
+              for (int u = 0; u < im->numCols; u++) {
+                if (isfinite(median->image->data.F32[v][u])) {
+
+                  im->data.F32[v][u] -= (median->image->data.F32[v][u]);//  - correction_median);
+                }
+              }
+            }
+            if (!psImageBackground(bgStats, NULL, ro->image, ro->mask, maskVal, rng)) {
+              psWarning("Unable to measure background for cell %d\n", j);
+              psErrorClear();
+              meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR;
+              continue;
+            }
+            psLogMsg("psModules.detrend",PS_LOG_DETAIL,"PATTERN.MEDIAN mode: %d Cell Statistics Post: %d %d %d %s %g %g\n",
+                     mode, j,x,y,cellName,bgStats->robustMedian,bgStats->robustStdev);
+            cellMedians->data.F32[j] = bgStats->robustMedian;
+            cellSigmas->data.F32[j]  = bgStats->robustStdev;
+          }
+        } // End subtraction block
+        // Write out the median?
+#if (1)
+        const char *expname = psMetadataLookupStr(NULL,chip->hdu->header,"FILENAME");
+        char filename[256];
+        pmCell *exCell2 = chip->cells->data[inputs->data.S32[0]];
+        snprintf(filename,sizeof(filename),
+                 "/home/panstarrs/watersc1/pattern_median.20130201/tests/%s.%d.%d.fits",
+                 expname,mode,i);
+        psFits *fits = psFitsOpen(filename,"w");
+        if (!fits) {
+          fprintf(stderr,"Couldn't write\n");
+        }
+        if (!psFitsWriteImage(fits,exCell2->hdu->header,median->image,0,NULL)) {
+          fprintf(stderr,"Couldn't write2\n");
+        }
+        if (!psFitsClose(fits)) {
+          fprintf(stderr,"Couldn't write3\n");
+        }
+#endif
+        
+        psFree(median->image);
+        psFree(median->mask);
+        psFree(median);
+      }
+    }
+    psFree(bgStats);
+    psFree(rng);
+    psFree(cellMedians);
+    psFree(cellSigmas);
+    psFree(meanMask);
+    psFree(stats);
+    psFree(inputs);
+    psFree(pixelData);
+    psFree(pixelMask);
+
+
+    return true;
+}