Changeset 42379 for trunk/psModules/src/detrend/pmPattern.c

Timestamp:

Feb 8, 2023, 11:54:08 AM (3 years ago)

Author:

eugene

Message:

merge from eam_branches/ipp-20220316. no_warn strncpy; avoid passing NULL to sprintf %s; add code for PATTERN_DEAD_CELLS; fix error in PSF residual image evaluation; drop detailed mask analysis for binned images (not actually used)

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/detrend/pmPattern.c (modified) (1 diff)

trunk/psModules/src/detrend/pmPattern.c

@@ -1005 +1005 @@
     psFree(meanMask);
 
+    return true;
+}
+
+// Compare the backs (cellBackgrounds) vector to each of the patterns loaded for this chip.
+// Choose the one that matches best & mask as appropriate
+// can we pass the backgrounds in using an image (8 x 8)
+bool pmPatternDeadCells(pmChip *chip, const psVector *backs, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_VECTOR_NON_NULL(backs, false);
+    PS_ASSERT_VECTOR_SIZE(backs, chip->cells->n, false);
+    PS_ASSERT_VECTOR_TYPE(backs, PS_TYPE_F32, false);
+
+    // Look for the dead cell pattern cube in the analysis metadata.
+    // NOTE: not all chips have the pattern, skip if not found.
+    bool mdok;
+    psImage *deadCellPattern = (psImage *) psMetadataLookupPtr (&mdok, chip->analysis, "PTN.DEAD.CELL");
+    if (!mdok) {
+        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No DEAD CELL pattern for chip, skipping\n");
+        return true;
+    }
+
+    int numCells = backs->n;            // Number of cells
+    int numColsD = deadCellPattern->numCols;
+    int numRowsD = deadCellPattern->numRows;
+
+    assert (backs->n == numRowsD);
+
+    // The background values need to be normalized (divide by the median).
+    // First extract the valid data values
+    psVector *valid = psVectorAllocEmpty(backs->n, PS_TYPE_F32); // Mean for each cell
+    for (int i = 0; i < backs->n; i++) {
+        float value = backs->data.F32[i];
+        if (!isfinite(value)) { backs->data.F32[i] = NAN; continue; }
+        if (value > 50000.0)  { backs->data.F32[i] = NAN; continue; } // XXX warning: hard-wired value
+        psVectorAppend (valid, value);
+    }
+    if (valid->n == 0) {
+        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No valid backgrounds, skipping\n");
+        psFree (valid);
+        return true;
+    }
+
+    // Second, calculate the median
+    psVectorSortInPlace (valid);
+    int midPt = valid->n / 2.0;
+    float median = valid->n % 2 ? valid->data.F32[midPt] : 0.5*(valid->data.F32[midPt] + valid->data.F32[midPt-1]);
+    psFree (valid);
+
+    // Finally, renormalize:
+    for (int i = 0; i < backs->n; i++) {
+        if (!isfinite(backs->data.F32[i])) { continue; }
+        backs->data.F32[i] /= median;
+    }
+
+    // there are 2 columns (value & mask) for each mode, plus the constant mode
+    int nModes = numColsD / 2 + 1;
+    
+    psVector *stdev = psVectorAlloc(nModes, PS_TYPE_F32); // Mean for each cell
+
+    // measure stdev for each comparison
+    for (int i = 0; i < nModes; i++) {
+        float Sum1 = 0.0;
+        float Sum2 = 0.0;
+        int   Npts = 0;
+
+        // choose the column with the background values for this mode
+        int nModeColumn = 2*(i - 1);
+
+        for (int j = 0; j < backs->n; j++) {
+
+            float valObs = backs->data.F32[j];
+            float valRef = (i == 0) ? 0.0 : deadCellPattern->data.F32[j][nModeColumn];
+
+            if (!isfinite(valObs)) continue;
+            if (!isfinite(valRef)) continue;
+
+            float dS = valObs - valRef;
+            Sum1 += dS;
+            Sum2 += dS*dS;
+            Npts ++;
+        }
+        if (Npts == 0) {
+            // is this is an error?
+            // no valid data, ignore this test
+            psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No valid backgrounds, skipping\n");
+            psFree (stdev);
+            return true;
+        }
+
+        float mean = Sum1 / Npts;
+        float sigma = sqrt(Sum2 / Npts - mean*mean);
+        stdev->data.F32[i] = sigma;
+        fprintf (stderr, "mode: %d, stdev: %f\n", i, sigma);
+    }
+
+    // loop over stdev and choose the lowest one
+    int   minI = 0;
+    float minV = stdev->data.F32[minI];
+
+    for (int i = 1; i < nModes; i++) {
+        if (stdev->data.F32[i] < minV) {
+            minI = i;
+            minV = stdev->data.F32[i];
+        }
+    }
+    psFree (stdev);
+
+    if (minI == 0) return true;
+
+    int nModeColumnMask = 2*(minI - 1) + 1;
+
+    // now mask bad cells
+    for (int i = 0; i < numCells; i++) {
+        if (deadCellPattern->data.F32[i][nModeColumnMask] > 0) continue;
+
+        pmCell *cell = chip->cells->data[i]; // Cell of interest
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+
+        psImage *mask = ro->mask; // mask of interest
+        int numCols = mask->numCols, numRows = mask->numRows; // Size of image
+        
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
+            }
+        }
+    }
     return true;
 }