Index: /trunk/psModules/src/detrend/pmShifts.c =================================================================== --- /trunk/psModules/src/detrend/pmShifts.c (revision 12010) +++ /trunk/psModules/src/detrend/pmShifts.c (revision 12010) @@ -0,0 +1,189 @@ +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#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 + + +// 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; + + +bool pmShiftsRead(psFits *fits, pmFPA *fpa, const pmFPAview *phuView) +{ + 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 + bool mdok; // Status of MD lookup + psMetadata *fileInfo = psMetadataLookupMetadata(&mdok, format, "FILE"); // FILE information from 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 + if (!mdok || !shiftsExt) { + // We read all the shifts that are present --- which is none. + return true; + } + + int numExt = psFitsGetSize(fits); // Number of extensions in 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"); + 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++; + } + 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); + } + psFree(namesIter); + psFree(names); + psFree(shifts); + + break; + } + } + + // reposition to the original position + return psFitsMoveExtNum(fits, origExt, false); +} Index: /trunk/psModules/src/detrend/pmShifts.h =================================================================== --- /trunk/psModules/src/detrend/pmShifts.h (revision 12010) +++ /trunk/psModules/src/detrend/pmShifts.h (revision 12010) @@ -0,0 +1,21 @@ +#ifndef PM_SHIFTS_H +#define PM_SHIFTS_H + +#include +#include "pmFPA.h" +#include "pmFPAview.h" + +/// Read any orthogonal transfer shifts extension in a FITS file +/// +/// 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 + ); + + + +#endif