Index: trunk/psModules/src/imcombine/Makefile.am =================================================================== --- trunk/psModules/src/imcombine/Makefile.am (revision 18140) +++ trunk/psModules/src/imcombine/Makefile.am (revision 18146) @@ -10,4 +10,5 @@ pmSubtraction.c \ pmSubtractionEquation.c \ + pmSubtractionIO.c \ pmSubtractionKernels.c \ pmSubtractionMask.c \ @@ -24,4 +25,5 @@ pmSubtraction.h \ pmSubtractionEquation.h \ + pmSubtractionIO.h \ pmSubtractionKernels.h \ pmSubtractionMask.h \ Index: trunk/psModules/src/imcombine/pmSubtractionIO.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 18146) +++ trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 18146) @@ -0,0 +1,284 @@ +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmSubtraction.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionMatch.h" + +#define ARRAY_BUFFER 16 // Number to add to array at a time + +// Names of FITS table columns +#define NAME_XMIN "XMIN" // Region of applicability: minimum x value +#define NAME_XMAX "XMAX" // Region of applicability: maximum x value +#define NAME_YMIN "YMIN" // Region of applicability: minimum y value +#define NAME_YMAX "YMAX" // Region of applicability: maximum y value +#define NAME_KERNEL "KERNEL" // Kernel description +#define NAME_TYPE "TYPE" // Kernel type +#define NAME_SIZE "SIZE" // Kernel half-size +#define NAME_INNER "INNER" // Size of inner region (only applicable for some kernel types) +#define NAME_SPATIAL "SPATIAL_ORDER" // Order of spatial polynomial +#define NAME_BG "BG_ORDER" // Order of background polynomial +#define NAME_MODE "MODE" // Matching mode +#define NAME_COLS "COLUMNS" // Number of columns +#define NAME_ROWS "ROWS" // Number of rows +#define NAME_SOL1 "SOLUTION_1" // Solution for convolving image 1 +#define NAME_SOL2 "SOLUTION_2" // Solution for convolving image 2 + + +// Generate an extension name for the subtraction kernel +static psString subtractionKernelExtname(const pmReadout *ro // Readout containing kernel + ) +{ + psAssert(ro, "Required parameter"); + + psString extname = NULL; // Name of FITS extension + pmCell *cell = ro->parent; // Parent cell + if (cell && cell->parent) { + pmChip *chip = cell->parent; // Parent chip + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + psAssert(chipName && cellName, "Chip and cell names should be defined in concepts."); + int roNum = -1; // Index of readout + for (int i = 0; i < cell->readouts->n && roNum == -1; i++) { + if (cell->readouts->data[i] == ro) { + roNum = i; + } + } + if (roNum != -1) { + psStringAppend(&extname, "%s_%s_%d", chipName, cellName, roNum); + } else { + psWarning("Unable to find readout number"); + } + } + + if (!extname) { + extname = psStringCopy("SUBTRACTION_KERNEL"); + } + + return extname; +} + +bool pmReadoutWriteSubtractionKernels(pmReadout *ro, psFits *fits) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PS_ASSERT_FITS_NON_NULL(fits, false); + + + // Extract the regions and solutions used in the image matching + psArray *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); + psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_REGION); + regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V); + } + psFree(iter); + } + if (regions->n == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No subtraction regions found."); + psFree(regions); + return false; + } + + psArray *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + // Set the normalisation dimensions, since these will be otherwise unavailable when reading the + // images by scans. + pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction + kernel->numCols = ro->image->numCols; + kernel->numRows = ro->image->numRows; + + kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); + } + psFree(iter); + } + + if (regions->n != kernels->n) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of regions (%ld) and kernels (%ld) don't match.\n", + regions->n, kernels->n); + psFree(regions); + psFree(kernels); + return false; + } + + // Format for writing a table + int num = regions->n; // Number of regions and kernels + psArray *rows = psArrayAlloc(num); // Array of FITS table rows + for (int i = 0; i < num; i++) { + psMetadata *row = psMetadataAlloc(); // Row of interest + rows->data[i] = psMemIncrRefCounter(row); + + psRegion *region = regions->data[i]; // Region of interest + psMetadataAddS32(row, PS_LIST_TAIL, NAME_XMIN, 0, "Applicability minimum x", region->x0); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_XMAX, 0, "Applicability maximum x", region->x1); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_YMIN, 0, "Applicability minimum y", region->y0); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_YMAX, 0, "Applicability maximum y", region->y1); + + pmSubtractionKernels *kernel = kernels->data[i]; // Kernel + psMetadataAddStr(row, PS_LIST_TAIL, NAME_KERNEL, 0, "Kernel description", kernel->description); + +#if 0 + psMetadataAddS32(row, PS_LIST_TAIL, NAME_TYPE, 0, "Kernel type (enum)", kernel->type); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_SIZE, 0, "Kernel half-size", kernel->size); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_INNER, 0, "Size of inner region", kernel->inner); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_SPATIAL, 0, "Polynomial order for spatial variations", + kernel->spatialOrder); +#endif + + psMetadataAddS32(row, PS_LIST_TAIL, NAME_BG, 0, "Polynomial order for background fitting", + kernel->bgOrder); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_MODE, 0, "Matching mode (enum)", kernel->mode); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_COLS, 0, "Number of columns", kernel->numCols); + psMetadataAddS32(row, PS_LIST_TAIL, NAME_ROWS, 0, "Number of rows", kernel->numRows); + if (kernel->mode == PM_SUBTRACTION_MODE_1 || kernel->mode == PM_SUBTRACTION_MODE_DUAL) { + psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL1, 0, "Solution vector 1", kernel->solution1); + } + if (kernel->mode == PM_SUBTRACTION_MODE_2 || kernel->mode == PM_SUBTRACTION_MODE_DUAL) { + psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL2, 0, "Solution vector 2", kernel->solution2); + } + } + psFree(regions); + psFree(kernels); + + psMetadata *header = psMetadataAlloc(); // Header for FITS file + + // CVS tags, used to identify the version of this file (in case incompatibilities are introduced) + psString cvsFile = psStringCopy("$RCSfile: pmSubtractionIO.c,v $"); + psString cvsRev = psStringCopy("$Revision: 1.1 $"); + psString cvsDate = psStringCopy("$Date: 2008-06-16 21:59:43 $"); + psStringSubstitute(&cvsFile, NULL, "RCSfile: "); + psStringSubstitute(&cvsRev, NULL, "Revision: "); + psStringSubstitute(&cvsDate, NULL, "Date: "); + + psString version = NULL; // Version information, for header + psStringAppend(&version, "%s %s %s", cvsFile, cvsRev, cvsDate); + psFree(cvsFile); + psFree(cvsRev); + psFree(cvsDate); + psStringSubstitute(&version, NULL, "$"); + psMetadataAddStr(header, PS_LIST_TAIL, "PSVERSION", 0, "S/W version", version); + psFree(version); + + psString extname = subtractionKernelExtname(ro); // Extension name + + if (!psFitsWriteTable(fits, header, rows, extname)) { + psError(PS_ERR_IO, false, "Unable to write subtraction kernel to FITS table."); + psFree(header); + psFree(rows); + psFree(extname); + return false; + } + + psFree(header); + psFree(rows); + psFree(extname); + + return true; +} + + +bool pmReadoutReadSubtractionKernels(pmReadout *ro, psFits *fits) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PS_ASSERT_FITS_NON_NULL(fits, false); + + psString extname = subtractionKernelExtname(ro); // Extension name + if (!psFitsMoveExtName(fits, extname)) { + psError(PS_ERR_IO, false, "Unable to move to subtraction kernel table."); + psFree(extname); + return false; + } + psFree(extname); + + psArray *table = psFitsReadTable(fits); // Table of interest + if (!table) { + psError(PS_ERR_IO, false, "Unable to read FITS table"); + return false; + } + + // Look up a column value for a row +#define TABLE_LOOKUP(TYPE, SUFFIX, TARGET, NAME) \ + TYPE TARGET; \ + { \ + bool mdok; \ + TARGET = psMetadataLookup##SUFFIX(&mdok, row, NAME); \ + if (!mdok) { \ + psError(PS_ERR_UNKNOWN, false, "Unable to find column %s in subtraction kernel table.", NAME); \ + psFree(table); \ + return false; \ + } \ + } + + for (int i = 0; i < table->n; i++) { + psMetadata *row = table->data[i]; // Table row + + TABLE_LOOKUP(int, S32, xMin, NAME_XMIN); + TABLE_LOOKUP(int, S32, xMax, NAME_XMAX); + TABLE_LOOKUP(int, S32, yMin, NAME_YMIN); + TABLE_LOOKUP(int, S32, yMax, NAME_YMAX); + + psRegion *region = psRegionAlloc(xMin, xMax, yMin, yMax); // Region of applicability + psMetadataAddPtr(ro->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_DATA_REGION, + "Subtraction region", region); + psFree(region); + + TABLE_LOOKUP(const char *, Str, description, NAME_KERNEL); + TABLE_LOOKUP(pmSubtractionMode, S32, mode, NAME_MODE); + +#if 0 + TABLE_LOOKUP(int, S32, size, NAME_SIZE); + TABLE_LOOKUP(int, S32, inner, NAME_INNER); + TABLE_LOOKUP(int, S32, spatial, NAME_SPATIAL); +#endif + + TABLE_LOOKUP(int, S32, bg, NAME_BG); + TABLE_LOOKUP(int, S32, numCols, NAME_COLS); + TABLE_LOOKUP(int, S32, numRows, NAME_ROWS); + + pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode); + kernels->numCols = numCols; + kernels->numRows = numRows; + + bool mdok; // Status of MD lookup + if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) { + kernels->solution1 = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, row, NAME_SOL1)); + if (!mdok) { + psError(PS_ERR_UNKNOWN, false, "Unable to find column %s in subtraction kernel table.", + NAME_SOL1); + psFree(kernels); + psFree(table); + return false; + } + } + if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) { + kernels->solution2 = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, row, NAME_SOL2)); + if (!mdok) { + psError(PS_ERR_UNKNOWN, false, "Unable to find column %s in subtraction kernel table.", + NAME_SOL2); + psFree(kernels); + psFree(table); + return false; + } + } + + psMetadataAddPtr(ro->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, PS_DATA_UNKNOWN, + "Subtraction kernels", kernels); + psFree(kernels); + } + psFree(table); + return true; +} Index: trunk/psModules/src/imcombine/pmSubtractionIO.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionIO.h (revision 18146) +++ trunk/psModules/src/imcombine/pmSubtractionIO.h (revision 18146) @@ -0,0 +1,22 @@ +#ifndef PM_SUBTRACTION_IO_H +#define PM_SUBTRACTION_IO_H + +#include + +#include +#include + +/// Write subtraction kernels within a readout to a FITS file +bool pmReadoutWriteSubtractionKernels( + pmReadout *readout, ///< Readout for which to write subtraction kernels (in analysis MD) + psFits *fits ///< FITS file to which to write + ); + + +bool pmReadoutReadSubtractionKernels( + pmReadout *readout, ///< Readout for which to read subtraction kernels (into analysis MD) + psFits *fits ///< FITS file to which to write + ); + + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 18140) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 18146) @@ -612,4 +612,95 @@ +// Intermediate string parsing functions required because of different APIs for strtol and strtof +static inline int parseStringInt(const char *string) +{ + return strtol(string, NULL, 10); +} +static inline float parseStringFloat(const char *string) +{ + return strtof(string, NULL); +} + + +// Parse a string of a number, up to some delimiter, and advance past the delimiter +#define PARSE_STRING_NUMBER(TARGET, STRING, DELIM, PARSEFUNC) { \ + char *start = STRING; /* Start of string */ \ + char *end = strchr(STRING, DELIM); /* End of string */ \ + if (!end) { \ + psAbort("End of string encountered"); \ + } \ + int stringSize = end - STRING; /* Size of string with value, NOT including \0 */ \ + char value[stringSize + 1]; /* String to parse */ \ + strncpy(value, start, stringSize); \ + value[stringSize] = '\0'; \ + TARGET = PARSEFUNC(value); \ + STRING += stringSize + 1; /* Advance past delimiter */ \ +} + + +pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder, + pmSubtractionMode mode) +{ + PS_ASSERT_STRING_NON_EMPTY(description, NULL); + + if (bgOrder != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Background order %d is not yet supported.", bgOrder); + return false; + } + + pmSubtractionKernelsType type = PM_SUBTRACTION_KERNEL_NONE; // Type of kernel + int size = 0; // Half-size of kernel + int spatialOrder = 0; // Order of spatial variations + const psVector *fwhms = NULL; // FWHM of Gaussians + const psVector *orders = NULL; // Polynomial order for each FWHM + int inner = 0; // Size of inner region + int binning = 0; // Binning to use + int ringsOrder = 0; // Polynomial order for rings + + if (strncmp(description, "ISIS", 4) == 0) { + // XXX Support for GUNK + if (strstr(description, "+POIS")) { + type = PM_SUBTRACTION_KERNEL_GUNK; + psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description); + } else { + type = PM_SUBTRACTION_KERNEL_ISIS; + char *ptr = (char*)description + 5; // Eat "ISIS(" + PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); + + // Count the number of Gaussians + int numGauss = 0; + for (char *string = ptr; string; string = strchr(string, '(')) { + numGauss++; + } + + fwhms = psVectorAlloc(numGauss, PS_TYPE_F32); + orders = psVectorAlloc(numGauss, PS_TYPE_S32); + + for (int i = 0; i < numGauss; i++) { + ptr++; // Eat the '(' + PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234," + PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)" + } + + ptr++; // Eat ',' + spatialOrder = parseStringInt(ptr); + } + } else if (strncmp(description, "RINGS", 5) == 0) { + type = PM_SUBTRACTION_KERNEL_RINGS; + char *ptr = (char*)description + 6; + PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); + PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt); + PARSE_STRING_NUMBER(ringsOrder, ptr, ',', parseStringInt); + PARSE_STRING_NUMBER(spatialOrder, ptr, ')', parseStringInt); + } else { + psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported."); + } + + + return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, + inner, binning, ringsOrder, mode); +} + + pmSubtractionKernelsType pmSubtractionKernelsTypeFromString(const char *type) { Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 18140) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 18146) @@ -177,4 +177,11 @@ ); +/// Generate a kernel using the description +pmSubtractionKernels *pmSubtractionKernelsFromDescription( + const char *description, ///< Description of kernel + int bgOrder, ///< Polynomial order for background fitting + pmSubtractionMode mode ///< Mode for subtraction + ); + /// Return the appropriate type from a string pmSubtractionKernelsType pmSubtractionKernelsTypeFromString(const char *string // String name for kernel type Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 18140) +++ trunk/psModules/src/psmodules.h (revision 18146) @@ -86,4 +86,5 @@ #include #include +#include #include #include