Changeset 16949 for trunk/psModules/src/detrend/pmDark.c

Timestamp:

Mar 12, 2008, 10:49:46 AM (18 years ago)

Author:

Paul Price

Message:

Allow the new-style darks to generate old-style-like darks. We normalise the dark images by a nominated concept (typically CELL.DARKTIME) and fit that. Allow *no* ordinates to be fitted in which case a constant is simply fit (as you usually would for a dark); this isn't the most efficient way of doing this, because it's solving the full least-squares equation instead of doing a straight average, but it's flexible. If there's a problem with speed, we can put in a special case later.

File:

• trunk/psModules/src/detrend/pmDark.c (modified) (27 diffs)

trunk/psModules/src/detrend/pmDark.c

@@ -1 +1 @@
 #include <stdio.h>
 #include <pslib.h>
+#include <string.h>
 
 #include "psPolynomialMD.h"
@@ -95 +96 @@
 
 
-bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, int iter, float rej,
-                   psMaskType maskVal)
+bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept,
+                   int iter, float rej, psMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(output, false);
@@ -109 +110 @@
     int numBadInputs = 0;               // Number of bad inputs
     psArray *values = psArrayAlloc(numInputs);
+    psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts
+    psVectorInit(roMask, 0);
+    psVector *norm = normConcept ? psVectorAlloc(numInputs, PS_TYPE_F32) : NULL; // Normalisations for each
     for (int i = 0; i < numInputs; i++) {
         values->data[i] = psVectorAlloc(numOrdinates, PS_TYPE_F32);
-    }
-    psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts
-    psVectorInit(roMask, 0);
+        if (norm) {
+            pmReadout *ro = inputs->data[i]; // Readout of interest
+            float normValue;            // Normalisation value
+            if (!ordinateLookup(&normValue, normConcept, false, NAN, NAN, ro)) {
+                psWarning("Unable to find value of %s for readout %d", normConcept, i);
+                roMask->data.U8[i] = 0xff;
+                norm->data.F32[i] = NAN;
+                numBadInputs++;
+                continue;
+            }
+            if (normValue == 0.0) {
+                psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i);
+                roMask->data.U8[i] = 0xff;
+                norm->data.F32[i] = NAN;
+                numBadInputs++;
+                continue;
+            }
+            norm->data.F32[i] = 1.0 / normValue;
+        }
+    }
     psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Orders for each concept
     for (int i = 0; i < numOrdinates; i++) {
@@ -135 +156 @@
 
         for (int j = 0; j < numInputs; j++) {
+            psVector *val = values->data[j]; // Value vector for readout
+            if (roMask->data.U8[j]) {
+                val->data.F32[i] = NAN;
+                continue;
+            }
+
             pmReadout *ro = inputs->data[j]; // Readout of interest
-            psVector *val = values->data[j]; // Value vector for readout
-
             float value = NAN;          // Value of ordinate
             if (!ordinateLookup(&value, ord->name, ord->scale, ord->min, ord->max, ro)) {
@@ -148 +173 @@
         }
     }
-    numInputs -= numBadInputs;
 
     if (psTraceGetLevel("psModules.detrend") > 9) {
@@ -162 +186 @@
     psFree(orders);
     int numTerms = poly->coeff->n;      // Number of terms in polynomial
-    if (numTerms > numInputs) {
+    if (numTerms > numInputs - numBadInputs) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%d) to fit polynomial terms (%d).",
-                numInputs, numTerms);
+                numInputs - numBadInputs, numTerms);
         psFree(values);
         psFree(roMask);
+        psFree(norm);
         return false;
     }
@@ -176 +201 @@
                                 inputs)) {
         psError(PS_ERR_UNKNOWN, false, "No valid input readouts.");
+        psFree(values);
+        psFree(roMask);
+        psFree(norm);
         return false;
     }
@@ -224 +252 @@
 
                 pixels->data.F32[r] = readout->image->data.F32[yIn][xIn];
+                if (norm) {
+                    pixels->data.F32[r] *= norm->data.F32[r];
+                }
                 if (readout->mask) {
                     mask->data.PS_TYPE_MASK_DATA[r] = readout->mask->data.PS_TYPE_MASK_DATA[yIn][xIn];
@@ -240 +271 @@
     }
 
+    psFree(norm);
     psFree(roMask);
     psFree(poly);
@@ -248 +280 @@
     psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                      PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
+    psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM,
+                     PS_META_REPLACE, "Dark normalisation", normConcept);
 
     for (int i = 0; i < numTerms; i++) {
@@ -286 +320 @@
         return false;
     }
+    bool mdok;                          // Status of MD lookup
+    psString normConcept = psMetadataLookupStr(&mdok, dark->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation
 
     int numOrdinates = ordinates->n;    // Number of ordinates
@@ -299 +335 @@
         values->data.F32[i] = value;
     }
+    float norm = NAN;                   // Normalisation value
+    if (normConcept) {
+        if (!ordinateLookup(&norm, normConcept, false, NAN, NAN, readout)) {
+            psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", normConcept);
+            psFree(values);
+            return false;
+        }
+    }
 
     psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Order for each polynomial
@@ -322 +366 @@
             }
             float value = psPolynomialMDEval(poly, values); // Value of dark current
+            if (normConcept) {
+                value *= norm;
+            }
             readout->image->data.F32[y][x] -= value;
             if (readout->mask && !isfinite(value)) {
@@ -352 +399 @@
 
     psArray *ordinates = NULL;      // Dark ordinates, to write
+    const char *normConcept = NULL;     // Normalisation concept
     psArray *chips = fpa->chips;    // Component chips
     for (int i = 0; i < chips->n; i++) {
@@ -359 +407 @@
             pmCell *cell = cells->data[j]; // Cell of interest
             bool mdok;              // Status of MD lookup
-            psArray *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES);
+            psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ords
             if (!mdok) {
                 continue;
             }
+            psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Norm
             if (ordinates) {
-                if (new != ordinates) {
+                if (newOrd != ordinates) {
                     psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
                     return false;
                 }
+                if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) ||
+                    (!normConcept && newNorm)) {
+                    psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells.");
+                    return false;
+                }
             } else {
-                ordinates = new;
+                ordinates = newOrd;
+                normConcept = newNorm;
             }
         }
@@ -379 +434 @@
     }
 
-    return pmDarkWrite(fits, hdu->header, ordinates);
+    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
 }
 
@@ -399 +454 @@
 
     psArray *ordinates = NULL;          // Dark ordinates, to write
+    const char *normConcept = NULL;     // Normalisation concept
     psArray *cells = chip->cells;       // Component cells
     for (int j = 0; j < cells->n; j++) {
         pmCell *cell = cells->data[j]; // Cell of interest
         bool mdok;                      // Status of MD lookup
-        psArray *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES);
+        psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ordinates
         if (!mdok) {
             continue;
         }
+        psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation
         if (ordinates) {
-            if (new != ordinates) {
+            if (newOrd != ordinates) {
                 psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
                 return false;
             }
+            if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) ||
+                (!normConcept && newNorm)) {
+                psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells.");
+                return false;
+            }
         } else {
-            ordinates = new;
+            ordinates = newOrd;
+            normConcept = newNorm;
         }
     }
@@ -422 +485 @@
     }
 
-    return pmDarkWrite(fits, hdu->header, ordinates);
+    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
 }
 
@@ -447 +510 @@
         return false;
     }
-
-    return pmDarkWrite(fits, hdu->header, ordinates);
-}
-
-bool pmDarkWrite(psFits *fits, const psMetadata *header, const psArray *ordinates)
+    bool mdok;                          // Status of MD lookup
+    psString normConcept = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation
+
+    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
+}
+
+bool pmDarkWrite(psFits *fits, psMetadata *header, const psArray *ordinates, const char *normConcept)
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
@@ -472 +537 @@
     }
 
-    // Format ordinates into FITS table
-    int numOrdinates = ordinates->n;// Number of ordinates
-    psArray *table = psArrayAlloc(numOrdinates); // FITS table, constructed from ordinates
-    for (int i = 0; i < ordinates->n; i++) {
-        pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
-        psMetadata *row = psMetadataAlloc(); // FITS table row
-        psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name);
-        psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order);
-        psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale);
-        psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MIN, 0, "Minimum value", ord->min);
-        psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MAX, 0, "Maximum value", ord->max);
-        table->data[i] = row;
-    }
-
-    if (!psFitsWriteTable(fits, header, table, PM_DARK_FITS_EXTNAME)) {
-        psError(PS_ERR_IO, false, "Unable to write dark ordinates.");
+    if (!psMemIncrRefCounter((psMetadata*)header)) {
+        header = psMetadataAlloc();
+    }
+    psMetadataAddStr(header, PS_LIST_TAIL, PM_DARK_HEADER_NORM, PS_META_REPLACE,
+                     "Dark normalisation concept", normConcept);
+
+    if (ordinates->n > 0) {
+        // Format ordinates into FITS table
+        int numOrdinates = ordinates->n;// Number of ordinates
+        psArray *table = psArrayAlloc(numOrdinates); // FITS table, constructed from ordinates
+        for (int i = 0; i < ordinates->n; i++) {
+            pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
+            psMetadata *row = psMetadataAlloc(); // FITS table row
+            psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name);
+            psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order);
+            psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale);
+            psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MIN, 0, "Minimum value", ord->min);
+            psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MAX, 0, "Maximum value", ord->max);
+            table->data[i] = row;
+        }
+
+        if (!psFitsWriteTable(fits, header, table, PM_DARK_FITS_EXTNAME)) {
+            psError(PS_ERR_IO, false, "Unable to write dark ordinates.");
+            psFree(table);
+            psFree(header);
+            return false;
+        }
         psFree(table);
-        return false;
-    }
-    psFree(table);
+    } else {
+        // No ordinates to write to a table, so write to a blank header.
+        if (!psFitsWriteBlank(fits, header, PM_DARK_FITS_EXTNAME)) {
+            psError(PS_ERR_IO, false, "Unable to write dark header.");
+            psFree(header);
+            return false;
+        }
+    }
+    psFree(header);
 
     return true;
@@ -507 +589 @@
     }
 
-    psArray *ordinates = pmDarkRead(fits);
+    psString normConcept = NULL;        // Normalisation concept
+    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
     if (!ordinates) {
         psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
@@ -521 +604 @@
             psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                              PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
+            psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
+                             "Dark normalisation", normConcept);
         }
     }
     psFree(ordinates);
+    psFree(normConcept);
 
     return true;
@@ -538 +624 @@
     }
 
-    psArray *ordinates = pmDarkRead(fits);
+    psString normConcept = NULL;        // Normalisation concept
+    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
     if (!ordinates) {
         psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
@@ -549 +636 @@
         psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                          PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
+        psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
+                         "Dark normalisation", normConcept);
     }
     psFree(ordinates);
+    psFree(normConcept);
 
     return true;
@@ -567 +657 @@
     }
 
-    psArray *ordinates = pmDarkRead(fits);
+    psString normConcept = NULL;        // Normalisation concept
+    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
     if (!ordinates) {
         psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
@@ -574 +665 @@
     psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                      PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
+    psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
+                     "Dark normalisation", normConcept);
     psFree(ordinates);
+    psFree(normConcept);
 
     return true;
 }
 
-psArray *pmDarkRead(psFits *fits)
-{
-    PS_ASSERT_FITS_NON_NULL(fits, false);
+psArray *pmDarkRead(psString *normConcept, psFits *fits)
+{
+    PS_ASSERT_PTR_NON_NULL(normConcept, NULL);
+    PS_ASSERT_PTR_NULL(*normConcept, NULL);
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
 
     if (!psFitsMoveExtName(fits, PM_DARK_FITS_EXTNAME)) {
@@ -589 +685 @@
     }
 
-    psArray *table = psFitsReadTable(fits); // FITS Table with ordinates
-    int numOrdinates = table->n;        // Number of ordinates
-    psArray *ordinates = psArrayAlloc(numOrdinates); // Parsed ordinates
-
-    for (int i = 0; i < numOrdinates; i++) {
-        psMetadata *row = table->data[i]; // Row of interest
-        const char *name = psMetadataLookupStr(NULL, row, PM_DARK_FITS_NAME); // Concept name
-        int order = psMetadataLookupS32(NULL, row, PM_DARK_FITS_ORDER); // Polynomial order
-        if (!name || order <= 0) {
-            psError(PS_ERR_UNKNOWN, false, "Bad value reading dark ordinates table.");
-            psFree(table);
-            psFree(ordinates);
-            return false;
-        }
-        pmDarkOrdinate *ord = pmDarkOrdinateAlloc(name, order); // Ordinate data
-        ord->scale = psMetadataLookupBool(NULL, row, PM_DARK_FITS_SCALE);
-        ord->min = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MIN);
-        ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX);
-
-        ordinates->data[i] = ord;
-    }
-    psFree(table);
+    psMetadata *header = psFitsReadHeader(NULL, fits); // Header
+    bool mdok;                          // Status of MD lookup
+    *normConcept = psMemIncrRefCounter(psMetadataLookupStr(&mdok, header, PM_DARK_HEADER_NORM));
+
+    psArray *ordinates = NULL;          // Dark ordinates to return
+
+    psFitsType type = psFitsGetExtType(fits); // Type of FITS extension
+    switch (type) {
+      case PS_FITS_TYPE_IMAGE: {
+          // Check that it's of zero size; otherwise we might have some conflict
+          int numCols = psMetadataLookupS32(&mdok, header, "NAXIS1");
+          int numRows = psMetadataLookupS32(&mdok, header, "NAXIS2");
+          if (numCols != 0 || numRows != 0) {
+              psError(PS_ERR_UNKNOWN, true, "Extension %s is not a DARK table.", PM_DARK_FITS_EXTNAME);
+              psFree(header);
+              return NULL;
+          }
+          // No ordinates fit --- only a constant term
+          ordinates = psArrayAlloc(0);
+          break;
+      }
+      case PS_FITS_TYPE_BINARY_TABLE:
+      case PS_FITS_TYPE_ASCII_TABLE: {
+          psArray *table = psFitsReadTable(fits); // FITS Table with ordinates
+          int numOrdinates = table->n;        // Number of ordinates
+          ordinates = psArrayAlloc(numOrdinates);
+
+          for (int i = 0; i < numOrdinates; i++) {
+              psMetadata *row = table->data[i]; // Row of interest
+              const char *name = psMetadataLookupStr(NULL, row, PM_DARK_FITS_NAME); // Concept name
+              int order = psMetadataLookupS32(NULL, row, PM_DARK_FITS_ORDER); // Polynomial order
+              if (!name || order <= 0) {
+                  psError(PS_ERR_UNKNOWN, false, "Bad value reading dark ordinates table.");
+                  psFree(table);
+                  psFree(ordinates);
+                  return false;
+              }
+              pmDarkOrdinate *ord = pmDarkOrdinateAlloc(name, order); // Ordinate data
+              ord->scale = psMetadataLookupBool(NULL, row, PM_DARK_FITS_SCALE);
+              ord->min = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MIN);
+              ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX);
+
+              ordinates->data[i] = ord;
+          }
+          psFree(table);
+          break;
+      }
+      default:
+        psError(PS_ERR_UNKNOWN, true, "Unrecognised FITS extension type.");
+        return NULL;
+    }
+    psFree(header);
 
     return ordinates;