Index: /trunk/psModules/src/detrend/pmShifts.c =================================================================== --- /trunk/psModules/src/detrend/pmShifts.c (revision 12190) +++ /trunk/psModules/src/detrend/pmShifts.c (revision 12191) @@ -12,8 +12,9 @@ #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,26 +23,33 @@ -// 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,7 +58,7 @@ // 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,26 +67,6 @@ } -// 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,5 +74,5 @@ 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,4 +81,5 @@ 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,5 +112,4 @@ 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,15 +117,66 @@ } + // 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,5 +193,5 @@ // 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,5 +203,5 @@ break; case PM_FPA_LEVEL_CELL: - singleShifts = otShiftsAlloc(); + singleShifts = pmShiftsAlloc(tRel, xyRel); break; default: @@ -177,7 +216,8 @@ 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,7 +228,8 @@ 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,21 +237,24 @@ } } - 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,5 +262,5 @@ } - 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,9 +271,9 @@ 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,9 +298,8 @@ 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,14 +336,14 @@ } 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,12 +361,13 @@ } 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,2 +384,97 @@ 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; +} Index: /trunk/psModules/src/detrend/pmShifts.h =================================================================== --- /trunk/psModules/src/detrend/pmShifts.h (revision 12190) +++ /trunk/psModules/src/detrend/pmShifts.h (revision 12191) @@ -4,4 +4,19 @@ #include #include "pmFPA.h" + +/// 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 + bool tRelative; ///< Are the time values relative (durations)? + bool xyRelative; ///< Are the shift (x,y) values relative to the previous position? +} pmShifts; + +/// Allocator for pmShifts +pmShifts *pmShiftsAlloc(bool tRel, ///< Are the time values relative (durations)? + bool xyRel ///< Are the shift (x,y) values relative to the previous position? + ); /// Read orthogonal transfer shifts table for a cell @@ -14,9 +29,20 @@ /// 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 +bool pmShiftsRead(const pmCell *cell, ///< Cell for which to search for shifts psFits *fits ///< FITS file in which to search for OT shifts extension ); +/// Generate a kernel for the cell from the orthogonal transfer shifts +/// +/// The kernel is saved in the analysis metadata +bool pmShiftsKernel(const pmCell *cell ///< Cell for which to generate kernel + ); + +/// Convolve a detrend with the appropriate orthogonal transfer convolution kernel from a science exposure. +/// +/// The kernel is generated (with pmShiftsKernel) if required. The image and mask are convolved with the +/// kernel (specified maskVal is smeared). The weight map is not convolved. +bool pmShiftsConvolve(pmReadout *detrend, const pmCell *source, psMaskType maskVal) #endif