Changeset 12098

Timestamp:

Feb 27, 2007, 4:14:42 PM (19 years ago)

Author:

Paul Price

Message:

Fixing up function to read and interpret orthogonal transfer shifts table. Compiles, but not yet tested.

Location:

trunk/psModules/src/detrend

Files:

• Makefile.am (modified) (2 diffs)
• pmShifts.c (modified) (4 diffs)
• pmShifts.h (modified) (1 diff)

trunk/psModules/src/detrend/Makefile.am

@@ -11 +11 @@
         pmBias.c \
         pmDetrendDB.c \
-        pmShutterCorrection.c
+        pmShutterCorrection.c \
+        pmShifts.c
 
 #       pmSkySubtract.c
@@ -23 +24 @@
         pmBias.h \
         pmDetrendDB.h \
-        pmShutterCorrection.h
+        pmShutterCorrection.h \
+        pmShifts.h
 
 #       pmSkySubtract.h

trunk/psModules/src/detrend/pmShifts.c

@@ -8 +8 @@
 #include <pslib.h>
 
+#include "pmFPALevel.h"
+#include "pmFPAUtils.h"
 #include "pmShifts.h"
 
@@ -13 +15 @@
 #define SHIFTS_BUFFER 100               // Buffer size for shifts
 #define ANALYSIS_NAME "OT.KERNEL"       // Name for kernel on the analysis metadata
+#define TRACE "psModules.detrend"       // Trace facility
+
+// XXX To do:
+// * Make the table column names configurable, by having a SHIFTS metadata in the camera format, with entries
+//   specifying the column names.
 
 
@@ -23 +30 @@
 } otShifts;
 
-
-bool pmShiftsRead(psFits *fits, pmFPA *fpa, const pmFPAview *phuView)
+// Allocator for otShifts
+static otShifts *otShiftsAlloc(void)
+{
+    otShifts *shifts = psAlloc(sizeof(otShifts));
+    shifts->x = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_S32);
+    shifts->y = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_S32);
+    shifts->t = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_F32);
+    shifts->num = 0;
+    return shifts;
+}
+
+// Look up the cell within a hash; supplement the hash with a new value if it doesn't exist
+static otShifts *cellVectors(psHash *shifts, // Hash of shifts
+                             const char *cellName // Key for hash
+    )
+{
+    assert(shifts);
+    assert(cellName);
+
+    // Find the appropriate cell
+    otShifts *vectors = psHashLookup(shifts, cellName);
+    if (!vectors) {
+        vectors = otShiftsAlloc();
+        psHashAdd(shifts, cellName, vectors);
+        psFree(vectors);            // Drop reference
+    }
+    return vectors;
+}
+
+// Generate a kernel and stuff it on the cell metadata
+static bool stuffKernel(pmCell *cell,     // Cell to which the shifts belong
+                        const otShifts *shifts // Shifts to apply
+    )
+{
+    assert(cell);
+    assert(shifts);
+
+    psKernel *kernel = psKernelGenerate(shifts->t, shifts->x, shifts->y, false, true); // Shift kernel
+    if (!kernel) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kernel from OT shifts");
+        return false;
+    }
+    psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, ANALYSIS_NAME, PS_DATA_KERNEL,
+                     "Orthogonal transfer kernel, calculated from shifts", kernel);
+    psFree(kernel);
+
+    return true;
+}
+
+
+bool pmShiftsRead(pmCell *cell, psFits *fits)
 {
     PS_ASSERT_PTR_NON_NULL(fits, false);
-    PS_ASSERT_PTR_NON_NULL(fpa, false);
-    PS_ASSERT_PTR_NON_NULL(phuView, false);
-
-    pmHDU *phu = pmFPAviewThisHDU(phuView, fpa); // The primary header
-    if (!phu) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Provided view does not correspond to an HDU.\n");
-        return false;
-    }
-
-    psMetadata *format = phu->format;   // Camera format
+    PS_ASSERT_PTR_NON_NULL(cell, false);
+
     bool mdok;                          // Status of MD lookup
-    psMetadata *fileInfo = psMetadataLookupMetadata(&mdok, format, "FILE"); // FILE information from format
+    psKernel *check = psMetadataLookupPtr(&mdok, cell->analysis, ANALYSIS_NAME); // Kernel, or NULL
+    if (check) {
+        psTrace(TRACE, 2, "Cell already has OT kernel present.\n");
+        return true;                    // No error
+    }
+    psTrace(TRACE, 2, "Reading FITS file for OT kernels.\n");
+
+    pmChip *chip = cell->parent;        // The parent chip
+    pmFPA *fpa = chip->parent;          // The parent FPA
+    pmHDU *phu = NULL;                  // The primary header
+    pmFPALevel phuLevel = PM_FPA_LEVEL_NONE; // Level of the PHU
+    long numChips = 0;                  // Number of chips below the PHU; for setting efficient hash size
+    long numCells = 0;                  // Number of cells below the PHU; for setting efficient hash size
+    if (fpa->hdu) {
+        phu = fpa->hdu;
+        // Count the cells
+        psArray *chips = fpa->chips;    // Array of chips
+        numChips = chips->n;
+        for (int i = 0; i < chips->n; i++) {
+            pmChip *chip = chips->data[i]; // Chip of interest
+            numCells += chip->cells->n;
+        }
+        numCells = (float)numCells / (float)numChips + 0.5; // Average number of cells per chip
+        phuLevel = PM_FPA_LEVEL_FPA;
+    }
+    if (!phu && chip->hdu) {
+        phu = chip->hdu;
+        numChips = 1;
+        numCells = chip->cells->n;
+        phuLevel = PM_FPA_LEVEL_CHIP;
+    }
+    if (!phu && cell->hdu) {
+        phu = cell->hdu;
+        numChips = 0;
+        numCells = 1;
+        phuLevel = PM_FPA_LEVEL_CELL;
+    }
+
+    if (!phu || phuLevel == PM_FPA_LEVEL_NONE) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Can't find the PHU.\n");
+        return false;
+    }
+
+    // Find out what to read
+    psMetadata *fileInfo = psMetadataLookupMetadata(&mdok, phu->format, "FILE"); // FILE in the format
     if (!mdok || !fileInfo) {
         psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find FILE information in camera format.\n");
@@ -49 +141 @@
     }
 
-    int numExt = psFitsGetSize(fits); // Number of extensions in FITS file
+    // Read the FITS file
     int origExt = psFitsGetExtNum(fits); // Original extension number; to preserve position
-
-    for (int i = 1; i < numExt; i++) {  // Note: shifts table cannot be in PHU!
-        psFitsMoveExtNum(fits, i, false);
-        const char *extname = psFitsGetExtName(fits); // Name of extension
-        if (strcmp(extname, shiftsExt) == 0) {
-            psArray *table = psFitsReadTable(fits); // The table of shifts
-            if (!table) {
-                psError(PS_ERR_IO, false, "Unable to read shifts table.\n");
+    if (!psFitsMoveExtName(fits, shiftsExt)) {
+        psError(PS_ERR_IO, false, "Unable to move to shifts extension %s", shiftsExt);
+        return false;
+    }
+    psArray *table = psFitsReadTable(fits); // The table of shifts
+    if (!table) {
+        psError(PS_ERR_IO, false, "Unable to read shifts table.\n");
+        return false;
+    }
+
+    // More sensible storage
+    otShifts *singleShifts = NULL;      // Shifts for a single cell
+    psHash *multipleShifts = NULL;      // Shifts for multiple cells, stored by cell name
+    switch (phuLevel) {
+      case PM_FPA_LEVEL_FPA:
+        multipleShifts = psHashAlloc(2 * numChips);
+        break;
+      case PM_FPA_LEVEL_CHIP:
+        multipleShifts = psHashAlloc(2 * numCells);
+        break;
+      case PM_FPA_LEVEL_CELL:
+        singleShifts = otShiftsAlloc();
+        break;
+      default:
+        psAbort("Should never get here.\n");
+    }
+
+    // Pull values out of the table into something a bit more sensible
+    for (int i = 0; i < table->n; i++) {
+        psMetadata *row = table->data[i]; // The row of interest
+
+        const char *chipName;           // Name of chip
+        if (phuLevel == PM_FPA_LEVEL_FPA) {
+            // Only care about the chip name if there's more than one chip
+            chipName = psMetadataLookupStr(&mdok, row, "Chip");
+            if (!mdok || !chipName) {
+                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Chip in row %d of shifts table\n", i);
+                psFree(multipleShifts);
+                psFree(table);
                 return false;
             }
-
-            // Pull values out of the table into something a bit more sensible
-            psHash *shifts = psHashAlloc(HASH_SIZE); // Shifts for each cell
-            for (int i = 0; i < table->n; i++) {
-                psMetadata *row = table->data[i]; // The row of interest
-                const char *name = psMetadataLookupStr(&mdok, row, "Cell"); // Name of cell
-                if (!mdok || !name) {
-                    psError(PS_ERR_UNEXPECTED_NULL, true,
-                            "Unable to find Cell in row %d of shifts table\n", i);
-                    psFree(shifts);
-                    psFree(table);
-                    return false;
-                }
-                float x = psMetadataLookupS32(&mdok, row, "X"); // Shift in x
-                if (!mdok) {
-                    psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find X in row %d of shifts table\n", i);
-                    psFree(shifts);
-                    psFree(table);
-                    return false;
-                }
-                float y = psMetadataLookupF32(&mdok, row, "Y"); // Shift in y
-                if (!mdok) {
-                    psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Y in row %d of shifts table\n", i);
-                    psFree(shifts);
-                    psFree(table);
-                    return false;
-                }
-                float t = psMetadataLookupF32(&mdok, row, "Time"); // Time of shift
-                if (!mdok) {
-                    psError(PS_ERR_UNEXPECTED_NULL, true,
-                            "Unable to find Time in row %d of shifts table\n", i);
-                    psFree(shifts);
-                    psFree(table);
-                    return false;
-                }
-
-                // Find the appropriate cell
-                otShifts *vectors = psHashLookup(shifts, name); // Shifts for the cell of interest
-                if (!vectors) {
-                    vectors = psAlloc(sizeof(otShifts));
-                    vectors->x = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_S32);
-                    vectors->y = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_S32);
-                    vectors->t = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_F32);
-                    vectors->num = 0;
-                    psHashAdd(shifts, name, vectors);
-                }
-
-                // Add the shift
-                psVectorExtend(cell->x, 1, SHIFTS_BUFFER);
-                psVectorExtend(cell->y, 1, SHIFTS_BUFFER);
-                psVectorExtend(cell->t, 1, SHIFTS_BUFFER);
-                vectors->x->data.S32[vectors->num] = (int)x;
-                vectors->y->data.S32[vectors->num] = (int)y;
-                vectors->t->data.F32[vectors->num] = t;
-                vectors->num++;
+        }
+        const char *cellName;           // Name of cell
+        if (phuLevel <= PM_FPA_LEVEL_CHIP) {
+            // Only care about the cell name if there's a chip
+            cellName = psMetadataLookupStr(&mdok, row, "Cell");
+            if (!mdok || !cellName) {
+                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Cell in row %d of shifts table\n", i);
+                psFree(multipleShifts);
+                psFree(table);
+                return false;
             }
+        }
+        float x = psMetadataLookupS32(&mdok, row, "X"); // Shift in x
+        if (!mdok) {
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find X in row %d of shifts table\n", i);
+            psFree(multipleShifts);
             psFree(table);
-
-            // Put the kernels into their own cells
-            psList *names = psHashKeyList(shifts); // List of cells
-            psListIterator *namesIter = psListIteratorAlloc(names, PS_LIST_HEAD, false); // Iterator for names
-            while ((name = psListGetAndIncrement(namesIter))) {
-                pmChip *chip = NULL;    // Chip of interest
-                pmCell *cell = NULL;    // Cell of interest
-                if (phuView->chip != -1) {
-                    chip = pmFPAviewThisChip(phuView, chip);
-                    cell = pmChipFindCell(chip, name);
-                } else {
-                    psList *chipCell = psStringSplit(name, ":|,", false); // The chip and cell
-                    if (chipCell->n != 2) {
-                        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                                "Cannot determine chip and cell names from %s\n", name);
-                        psFree(chipCell);
-                        psFree(namesIter);
-                        psFree(names);
-                        psFree(shifts);
-                        return false;
-                    }
-                    const char *chipName = psListGet(chipCell, PS_LIST_HEAD); // Name of chip
-                    const char *cellName = psListGet(chipCell, PS_LIST_TAIL); // Name of cell
-                    chip = pmFPAFindChip(fpa, chipName);
-                    if (!chip) {
-                        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find chip %s", chipName);
-                        psFree(chipCell);
-                        psFree(namesIter);
-                        psFree(names);
-                        psFree(shifts);
-                        return false;
-                    }
-                    cell = pmChipFindCell(chip, cellName);
-                    psFree(chipCell);
-                }
-
-                if (!cell) {
-                    psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find chip %s", chipName);
-                    psFree(namesIter);
-                    psFree(names);
-                    psFree(shifts);
-                    return false;
-                }
-
-                otShifts *vectors = psHashLookup(shifts, name); // Shifts for the cell of interest
-                psKernel *kernel = psKernelGenerate(vectors->t, vectors->x, vectors->y,
-                                                    false, true); // Shift kernel
-                if (!kernel) {
-                    psError(PS_ERR_UNEXPECTED_NULL, false,
-                            "Unable to generate kernel from OT shifts for %s", name);
-                    psFree(shifts);
-                    psFree(namesIter);
-                    psFree(names);
-                    return false;
-                }
-                psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, ANALYSIS_NAME, PS_DATA_KERNEL,
-                                 "Orthogonal transfer kernel, calculated from shifts", kernel);
-                psFree(kernel);
+            return false;
+        }
+        float y = psMetadataLookupF32(&mdok, row, "Y"); // Shift in y
+        if (!mdok) {
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Y in row %d of shifts table\n", i);
+            psFree(multipleShifts);
+            psFree(table);
+            return false;
+        }
+        float t = psMetadataLookupF32(&mdok, row, "Time"); // Time of shift
+        if (!mdok) {
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Time in row %d of shifts table\n", i);
+            psFree(multipleShifts);
+            psFree(table);
+            return false;
+        }
+
+        otShifts *shifts = NULL;        // Shifts for the cell of interest
+        switch (phuLevel) {
+          case PM_FPA_LEVEL_FPA: {
+              psHash *cells = psHashLookup(multipleShifts, chipName); // Hash of cells
+              if (!cells) {
+                  cells = psHashAlloc(numCells);
+                  psHashAdd(multipleShifts, chipName, cells);
+                  psFree(cells);          // Drop reference
+              }
+              shifts = cellVectors(cells, cellName);
+              break;
+          }
+          case PM_FPA_LEVEL_CHIP:
+            shifts = cellVectors(multipleShifts, cellName);
+            break;
+          case PM_FPA_LEVEL_CELL:
+            shifts = singleShifts;
+            break;
+          default:
+            psAbort("Should never get here.\n");
+        }
+
+        // Add the shift
+        psVectorExtend(shifts->x, 1, SHIFTS_BUFFER);
+        psVectorExtend(shifts->y, 1, SHIFTS_BUFFER);
+        psVectorExtend(shifts->t, 1, SHIFTS_BUFFER);
+        shifts->x->data.S32[shifts->num] = (int)x;
+        shifts->y->data.S32[shifts->num] = (int)y;
+        shifts->t->data.F32[shifts->num] = t;
+        shifts->num++;
+    }
+    psFree(table);
+
+    // Put the kernels into their own cells
+    if (phuLevel == PM_FPA_LEVEL_CELL) {
+        // Only a single cell
+        if (!stuffKernel(cell, singleShifts)) {
+            psFree(singleShifts);
+            return false;
+        }
+        psFree(singleShifts);
+    } else {
+        psList *names = psHashKeyList(multipleShifts); // List of hash keys (chip/cell names)
+        psListIterator *namesIter = psListIteratorAlloc(names, PS_LIST_HEAD, false); // Iterator for names
+        const char *name;               // Name, from iteration
+        while ((name = psListGetAndIncrement(namesIter))) {
+            switch (phuLevel) {
+              case PM_FPA_LEVEL_FPA: {
+                  int chipNum = pmFPAFindChip(fpa, name); // Number of chip of interest
+                  if (chipNum < 0) {
+                      psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find chip %s", name);
+                      psFree(namesIter);
+                      psFree(names);
+                      psFree(multipleShifts);
+                      return false;
+                  }
+                  pmChip *chip = fpa->chips->data[chipNum]; // Chip of interest
+
+                  // Loop over component cells
+                  psHash *cells = psHashLookup(multipleShifts, name); // Hash of cells
+                  psList *cellNames = psHashKeyList(multipleShifts); // List of hash keys (cell names)
+                  psListIterator *cellNamesIter = psListIteratorAlloc(cellNames, PS_LIST_HEAD, false);
+                  const char *cellName;       // Cell name, from iteration
+                  while ((cellName = psListGetAndIncrement(cellNamesIter))) {
+                      int cellNum = pmChipFindCell(chip, cellName); // Number of cell of interest
+                      if (!cell) {
+                          psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find cell %s", cellName);
+                          psFree(cellNamesIter);
+                          psFree(cellNames);
+                          psFree(namesIter);
+                          psFree(names);
+                          psFree(multipleShifts);
+                          return false;
+                      }
+                      pmCell *cell = chip->cells->data[cellNum]; // Cell of interest
+
+                      otShifts *vectors = psHashLookup(cells, cellName); // Shifts for the cell of interest
+                      if (!stuffKernel(cell, vectors)) {
+                          psFree(cellNamesIter);
+                          psFree(cellNames);
+                          psFree(namesIter);
+                          psFree(names);
+                          psFree(multipleShifts);
+                          return false;
+                      }
+                  }
+                  psFree(cellNamesIter);
+                  psFree(cellNames);
+                  break;
+              }
+              case PM_FPA_LEVEL_CHIP: {
+                  int cellNum = pmChipFindCell(chip, name); // Number of cell of interest
+                  if (!cell) {
+                      psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find cell %s", name);
+                      psFree(namesIter);
+                      psFree(names);
+                      psFree(multipleShifts);
+                      return false;
+                  }
+                  pmCell *cell = chip->cells->data[cellNum]; // Cell of interest
+
+                  otShifts *vectors = psHashLookup(multipleShifts, name); // Shifts for this cell
+                  if (!stuffKernel(cell, vectors)) {
+                      psFree(namesIter);
+                      psFree(names);
+                      psFree(multipleShifts);
+                      return false;
+                  }
+                  break;
+              }
+              default:
+                psAbort("Should never get here.\n");
             }
-            psFree(namesIter);
-            psFree(names);
-            psFree(shifts);
-
-            break;
-        }
-    }
-
-    // reposition to the original position
+        }
+        psFree(namesIter);
+        psFree(names);
+        psFree(multipleShifts);
+    }
+
+    // Go back to the original position in the FITS file
     return psFitsMoveExtNum(fits, origExt, false);
 }

trunk/psModules/src/detrend/pmShifts.h

@@ -4 +4 @@
 #include <pslib.h>
 #include "pmFPA.h"
-#include "pmFPAview.h"
 
-/// Read any orthogonal transfer shifts extension in a FITS file
+/// Read orthogonal transfer shifts table for a cell
 ///
-/// Given a FITS file (typically newly opened ready for processing), this function searches for the extension
-/// specified by the SHIFTS keyword within the FILE information in the camera format.  If the extension is
-/// found, the shifts table is read and translated into kernels which are placed in the analysis metadata of
-/// each cell.
-bool pmShiftsRead(psFits *fits,         ///< FITS file in which to search for OT shifts extension
-                  pmFPA *fpa,           ///< FPA containing the cells that own the shifts
-                  const pmFPAview *phuView ///< View of the PHU; specifies which chip(s) owns the cells
-    );
+/// Given a cell, this function searches for the orthogonal transfer shifts for this cell in the supplied FITS
+/// file.  The FITS extension containing the shifts table is specified by the SHIFTS keyword within the FILE
+/// information in the camera format.  If the extension is found, the shifts table is read and translated into
+/// kernels which are placed in the analysis metadata of each cell.  Note that this operation is performed on
+/// all cells within the FITS file (or at least, those contained within the shifts table), not just the cell
+/// provided --- if we have to read the whole table, we may as well translate the whole lot.  If a kernel is
+/// already present in the cell analysis metadata, the function returns true without doing any work.
+bool pmShiftsRead(pmCell *cell,         ///< Cell for which to search for shifts
+                  psFits *fits          ///< FITS file in which to search for OT shifts extension
+                  );