Index: trunk/psModules/src/detrend/pmShifts.c =================================================================== --- trunk/psModules/src/detrend/pmShifts.c (revision 12010) +++ trunk/psModules/src/detrend/pmShifts.c (revision 12098) @@ -8,4 +8,6 @@ #include +#include "pmFPALevel.h" +#include "pmFPAUtils.h" #include "pmShifts.h" @@ -13,4 +15,9 @@ #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,20 +30,105 @@ } 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 +141,201 @@ } - 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); }