Changeset 12191 for trunk/psModules/src/detrend/pmShifts.c

Timestamp:

Mar 2, 2007, 2:33:50 PM (19 years ago)

Author:

Paul Price

Message:

Splitting out reading shifts table and constructing the kernel (we might want access to the shifts table for other purposes besides convolution, like constructing a single image from the video data). Making shifts type public (pmShifts). Adding function to do the convolution.

File:

• trunk/psModules/src/detrend/pmShifts.c (modified) (19 diffs)

trunk/psModules/src/detrend/pmShifts.c

@@ -12 +12 @@
 #include "pmShifts.h"
 
-#define HASH_SIZE 100                   // Size of hash with cells; should be > numCells for efficiency
 #define SHIFTS_BUFFER 100               // Buffer size for shifts
-#define ANALYSIS_NAME "OT.KERNEL"       // Name for kernel on the analysis metadata
+#define TABLE_NAME "SHIFTS.TABLE"       // Name for table on the analysis metadata
+#define KERNEL_NAME "SHIFTS.KERNEL"     // Name for kernel on the analysis metadata
 #define TRACE "psModules.detrend"       // Trace facility
+#define FFT_SIZE 25                     // Size at which we use FFT instead of direct convolution
 
 // XXX To do:
@@ -22 +23 @@
 
 
-// Shifts due to orthogonal transfer
-typedef struct {
-    psVector *x;                        // Shifts in x
-    psVector *y;                        // Shifts in y
-    psVector *t;                        // Times of shifts
-    long num;                           // Number of values
-} otShifts;
-
-// Allocator for otShifts
-static otShifts *otShiftsAlloc(void)
-{
-    otShifts *shifts = psAlloc(sizeof(otShifts));
+static void shiftsFree(pmShifts *shifts)
+{
+    psFree(shifts->x);
+    psFree(shifts->y);
+    psFree(shifts->t);
+    return;
+}
+
+pmShifts *pmShiftsAlloc(bool tRel, bool xyRel)
+{
+    pmShifts *shifts = psAlloc(sizeof(pmShifts));
+    psMemSetDeallocator(shifts, (psFreeFunc)shiftsFree);
+
     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;
+
+    // Suitable defaults
+    shifts->tRelative = tRel;
+    shifts->xyRelative = xyRel;
+
     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
+static pmShifts *cellVectors(psHash *shifts, // Hash of shifts
+                             const char *cellName, // Key for hash
+                             bool tRel, bool xyRel // Are the shifts relative?
     )
 {
@@ -50 +58 @@
 
     // Find the appropriate cell
-    otShifts *vectors = psHashLookup(shifts, cellName);
+    pmShifts *vectors = psHashLookup(shifts, cellName);
     if (!vectors) {
-        vectors = otShiftsAlloc();
+        vectors = pmShiftsAlloc(tRel, xyRel);
         psHashAdd(shifts, cellName, vectors);
         psFree(vectors);            // Drop reference
@@ -59 +67 @@
 }
 
-// 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)
+
+bool pmShiftsRead(const pmCell *cell, psFits *fits)
 {
     PS_ASSERT_PTR_NON_NULL(fits, false);
@@ -86 +74 @@
 
     bool mdok;                          // Status of MD lookup
-    psKernel *check = psMetadataLookupPtr(&mdok, cell->analysis, ANALYSIS_NAME); // Kernel, or NULL
+    pmShifts *check = psMetadataLookupPtr(&mdok, cell->analysis, TABLE_NAME); // Table, or NULL
     if (check) {
         psTrace(TRACE, 2, "Cell already has OT kernel present.\n");
@@ -93 +81 @@
     psTrace(TRACE, 2, "Reading FITS file for OT kernels.\n");
 
+    // Determine camera layout
     pmChip *chip = cell->parent;        // The parent chip
     pmFPA *fpa = chip->parent;          // The parent FPA
@@ -123 +112 @@
         phuLevel = PM_FPA_LEVEL_CELL;
     }
-
     if (!phu || phuLevel == PM_FPA_LEVEL_NONE) {
         psError(PS_ERR_UNEXPECTED_NULL, true, "Can't find the PHU.\n");
@@ -129 +117 @@
     }
 
+
     // 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");
-        return false;
-    }
-    const char *shiftsExt = psMetadataLookupStr(&mdok, fileInfo, "SHIFTS"); // Name of shifts extension
+    psMetadata *shiftsInfo = psMetadataLookupMetadata(&mdok, phu->format, "SHIFTS"); // Shifts information
+    if (!mdok || !shiftsInfo) {
+        // We have read all the shifts that we have been told about --- which is none.
+        return true;
+    }
+    const char *shiftsExt = psMetadataLookupStr(&mdok, shiftsInfo, "EXTENSION"); // Extension name for shifts
     if (!mdok || !shiftsExt) {
-        // We read all the shifts that are present --- which is none.
-        return true;
-    }
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Can't find EXTENSION name in SHIFTS information from camera format.");
+        return false;
+    }
+    const char *chipCol = NULL;         // Column name for chip
+    if (phuLevel == PM_FPA_LEVEL_FPA) {
+        chipCol = psMetadataLookupStr(&mdok, shiftsInfo, "CHIP");
+        if (!mdok || !chipCol) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Can't find CHIP column name in SHIFTS information from camera format.");
+            return false;
+        }
+    }
+    const char *cellCol = NULL;         // Column name for cell
+    if (phuLevel <= PM_FPA_LEVEL_CHIP) {
+        cellCol = psMetadataLookupStr(&mdok, shiftsInfo, "CELL");
+        if (!mdok || !chipCol) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Can't find CELL column name in SHIFTS information from camera format.");
+            return false;
+        }
+    }
+    const char *tCol = psMetadataLookupStr(&mdok, shiftsInfo, "T");
+    if (!mdok || !tCol) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Can't find T column name in SHIFTS information from camera format.");
+        return false;
+    }
+    const char *xCol = psMetadataLookupStr(&mdok, shiftsInfo, "X");
+    if (!mdok || !xCol) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Can't find X column name in SHIFTS information from camera format.");
+        return false;
+    }
+    const char *yCol = psMetadataLookupStr(&mdok, shiftsInfo, "Y");
+    if (!mdok || !yCol) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Can't find Y column name in SHIFTS information from camera format.");
+        return false;
+    }
+
+    bool tRel = psMetadataLookupBool(&mdok, shiftsInfo, "TRELATIVE");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Can't find TRELATIVE in SHIFTS information from camera format.");
+        return false;
+    };
+    bool xyRel = psMetadataLookupStr(&mdok, shiftsInfo, "XYRELATIVE");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Can't find XYRELATIVE in SHIFTS information from camera format.");
+        return false;
+    };
 
     // Read the FITS file
@@ -154 +193 @@
 
     // More sensible storage
-    otShifts *singleShifts = NULL;      // Shifts for a single cell
+    pmShifts *singleShifts = NULL;      // Shifts for a single cell
     psHash *multipleShifts = NULL;      // Shifts for multiple cells, stored by cell name
     switch (phuLevel) {
@@ -164 +203 @@
         break;
       case PM_FPA_LEVEL_CELL:
-        singleShifts = otShiftsAlloc();
+        singleShifts = pmShiftsAlloc(tRel, xyRel);
         break;
       default:
@@ -177 +216 @@
         if (phuLevel == PM_FPA_LEVEL_FPA) {
             // Only care about the chip name if there's more than one chip
-            chipName = psMetadataLookupStr(&mdok, row, "Chip");
+            chipName = psMetadataLookupStr(&mdok, row, chipCol);
             if (!mdok || !chipName) {
-                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Chip in row %d of shifts table\n", i);
+                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
+                        chipCol, i);
                 psFree(multipleShifts);
                 psFree(table);
@@ -188 +228 @@
         if (phuLevel <= PM_FPA_LEVEL_CHIP) {
             // Only care about the cell name if there's a chip
-            cellName = psMetadataLookupStr(&mdok, row, "Cell");
+            cellName = psMetadataLookupStr(&mdok, row, cellCol);
             if (!mdok || !cellName) {
-                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Cell in row %d of shifts table\n", i);
+                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
+                        cellCol, i);
                 psFree(multipleShifts);
                 psFree(table);
@@ -196 +237 @@
             }
         }
-        float x = psMetadataLookupS32(&mdok, row, "X"); // Shift in x
+        float x = psMetadataLookupS32(&mdok, row, xCol); // Shift in x
         if (!mdok) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find X in row %d of shifts table\n", i);
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
+                    xCol, i);
             psFree(multipleShifts);
             psFree(table);
             return false;
         }
-        float y = psMetadataLookupF32(&mdok, row, "Y"); // Shift in y
+        float y = psMetadataLookupF32(&mdok, row, yCol); // Shift in y
         if (!mdok) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Y in row %d of shifts table\n", i);
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
+                    yCol, i);
             psFree(multipleShifts);
             psFree(table);
             return false;
         }
-        float t = psMetadataLookupF32(&mdok, row, "Time"); // Time of shift
+        float t = psMetadataLookupF32(&mdok, row, tCol); // Time of shift
         if (!mdok) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Time in row %d of shifts table\n", i);
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
+                    tCol, i);
             psFree(multipleShifts);
             psFree(table);
@@ -218 +262 @@
         }
 
-        otShifts *shifts = NULL;        // Shifts for the cell of interest
+        pmShifts *shifts = NULL;        // Shifts for the cell of interest
         switch (phuLevel) {
           case PM_FPA_LEVEL_FPA: {
@@ -227 +271 @@
                   psFree(cells);          // Drop reference
               }
-              shifts = cellVectors(cells, cellName);
+              shifts = cellVectors(cells, cellName, tRel, xyRel);
               break;
           }
           case PM_FPA_LEVEL_CHIP:
-            shifts = cellVectors(multipleShifts, cellName);
+            shifts = cellVectors(multipleShifts, cellName, tRel, xyRel);
             break;
           case PM_FPA_LEVEL_CELL:
@@ -254 +298 @@
     if (phuLevel == PM_FPA_LEVEL_CELL) {
         // Only a single cell
-        if (!stuffKernel(cell, singleShifts)) {
-            psFree(singleShifts);
-            return false;
-        }
+        psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, TABLE_NAME, PS_DATA_KERNEL,
+                         "Orthogonal transfer shifts", singleShifts);
         psFree(singleShifts);
+        return true;
     } else {
         psList *names = psHashKeyList(multipleShifts); // List of hash keys (chip/cell names)
@@ -293 +336 @@
                       }
                       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;
+                      if (psMetadataLookup(cell->analysis, TABLE_NAME)) {
+                          // Already has a shifts table, for some reason
+                          psWarning("Chip %s, cell %s already has a shifts table --- overwriting\n",
+                                    name, cellName);
                       }
+
+                      pmShifts *vectors = psHashLookup(cells, cellName); // Shifts for the cell of interest
+                      psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, TABLE_NAME,
+                                       PS_DATA_KERNEL | PS_META_REPLACE,
+                                       "Orthogonal transfer shifts", vectors);
                   }
                   psFree(cellNamesIter);
@@ -318 +361 @@
                   }
                   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;
+                  if (psMetadataLookup(cell->analysis, TABLE_NAME)) {
+                      // Already has a shifts table, for some reason
+                      psWarning("Cell %s already has a shifts table --- overwriting\n", name);
                   }
+
+                  pmShifts *vectors = psHashLookup(multipleShifts, name); // Shifts for this cell
+                  psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, TABLE_NAME,
+                                   PS_DATA_KERNEL | PS_META_REPLACE,
+                                   "Orthogonal transfer shifts", vectors);
                   break;
               }
@@ -340 +384 @@
     return psFitsMoveExtNum(fits, origExt, false);
 }
+
+
+// Generate a kernel and stuff it on the cell metadata
+bool pmShiftsKernel(const pmCell *cell     // Cell to which the shifts belong
+    )
+{
+    PS_ASSERT_PTR(cell, false);
+
+    bool mdok;                          // Status of MD lookup
+    pmShifts *shifts = psMetadataLookupPtr(&mdok, cell->analysis, TABLE_NAME);
+    if (!shifts) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find shifts table for cell.\n");
+        return false;
+    }
+
+    psKernel *kernel = psKernelGenerate(shifts->t, shifts->x, shifts->y,
+                                        shifts->tRelative, shifts->xyRelative); // 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, KERNEL_NAME, PS_DATA_KERNEL,
+                     "Orthogonal transfer kernel, calculated from shifts", kernel);
+    psFree(kernel);
+
+    return true;
+}
+
+bool pmShiftsConvolve(pmReadout *detrend, const pmCell *source, psMaskType maskVal)
+{
+    PS_ASSERT_PTR(detrend, false);
+    PS_ASSERT_PTR(source, false);
+
+    bool mdok;                          // Status of MD lookup
+    psKernel *kernel = psMetadataLookupPtr(&mdok, source->analysis, KERNEL_NAME);
+    if (!kernel) {
+        // Maybe they just forgot to generate the kernel with pmShiftsKernel
+        if (psMetadataLookup(source->analysis, TABLE_NAME)) {
+            if (!pmShiftsKernel(source)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to generate shifts kernel.");
+                return false;
+            }
+            // Hopefully it's there now
+            kernel = psMetadataLookupPtr(&mdok, source->analysis, KERNEL_NAME);
+            if (!kernel) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to find shifts kernel.");
+                return false;
+            }
+        } else {
+            psError(PS_ERR_UNKNOWN, false, "Unable to find shifts kernel or shifts table.");
+            return false;
+        }
+    }
+
+    if (detrend->image) {
+#if 1
+        // Always do direct convolution (no fuss with edge effects)
+        psImage *convolved = psImageConvolveDirect(detrend->image, kernel);
+#else
+        // Kernel size-dependent convolution --- if it's big, use the FFT
+        int xSize = kernel->xMax - kernel->xMin; // Kernel size in x
+        int ySize = kernel->yMax - kernel->yMin; // Kernel size in y
+        psImage *convolved;
+        if (xSize * ySize < FFT_SIZE * FFT_SIZE) {
+            convolved = psImageConvolveDirect(detrend->image, kernel);
+        } else {
+            // This is a little dodgy --- making choices about parameters without the user's input
+            psStats *stats = psImageStats(PS_STAT_ROBUST_MEDIAN);
+            stats->nSubsample = 10000;
+            psImageBackground(stats, detrend->image, detrend->mask, maskVal, NULL);
+            convolved = psImageConvolveFFT(detrend->image, kernel, stats->robustMedian);
+        }
+#endif
+        if (!convolved) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve detrend image.");
+            return false;
+        }
+        psFree(detrend->image);
+        detrend->image = convolved;
+    }
+
+    // Purposely ignoring the weight map --- don't care about the weight map for a detrend image
+
+    if (maskVal && detrend->mask) {
+        psImage *convolved = psImageConvolveMask(detrend->mask, maskVal, kernel);
+        if (!convolved) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve detrend mask.");
+            return false;
+        }
+        psFree(detrend->mask);
+        detrend->mask = convolved;
+    }
+
+    return true;
+}