Index: branches/eam_branches/psModules.stack.20100120/src/camera/pmFPAMaskWeight.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/camera/pmFPAMaskWeight.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/camera/pmFPAMaskWeight.c (revision 26747) @@ -370,8 +370,34 @@ psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts - int numCols = image->numCols, numRows = image->numRows; // Size of image - int numPix = numCols * numRows; // Number of pixels - int num = PS_MAX(sample, numPix); // Number we care about + + int xMin, xMax, yMin, yMax; // Bounds of image + if (mask) { + xMin = numCols; + xMax = 0; + yMin = numRows; + yMax = 0; + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) { + continue; + } + xMin = PS_MIN(xMin, x); + xMax = PS_MAX(xMax, x); + yMin = PS_MIN(yMin, y); + yMax = PS_MAX(yMax, y); + } + } + } else { + xMin = 0; + xMax = numCols; + yMin = 0; + yMax = numRows; + } + + int xNum = xMax - xMin, yNum = yMax - yMin; // Number of pixels + + int numPix = xNum * yNum; // Number of pixels + int num = PS_MAX(sample, numPix); // Number we care about psVector *signoise = psVectorAllocEmpty(num, PS_TYPE_F32); // Signal-to-noise values @@ -379,6 +405,6 @@ // We have an image smaller than Nsubset, so just loop over the image pixels int index = 0; // Index for vector - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { + for (int y = yMin; y < yMax; y++) { + for (int x = xMin; x < xMax; x++) { if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) { @@ -397,6 +423,6 @@ // Pixel coordinates int pixel = numPix * psRandomUniform(rng); - int x = pixel % numCols; - int y = pixel / numCols; + int x = xMin + pixel % xNum; + int y = yMin + pixel / xNum; if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || @@ -433,5 +459,7 @@ } - psBinaryOp(variance, variance, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32)); + psImage *subImage = psImageSubset(variance, psRegionSet(xMin, xMax, yMin, yMax)); // Smaller image + psBinaryOp(subImage, subImage, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32)); + psFree(subImage); pmHDU *hdu = pmHDUFromReadout(readout); // HDU for readout Index: branches/eam_branches/psModules.stack.20100120/src/config/Makefile.am =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/config/Makefile.am (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/config/Makefile.am (revision 26747) @@ -32,4 +32,5 @@ pmConfigDump.c \ pmConfigRun.c \ + pmConfigRecipeValue.c \ pmVersion.c \ pmErrorCodes.c @@ -43,4 +44,5 @@ pmConfigDump.h \ pmConfigRun.h \ + pmConfigRecipeValue.h \ pmVersion.h \ pmErrorCodes.h Index: branches/eam_branches/psModules.stack.20100120/src/config/pmConfigRecipeValue.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/config/pmConfigRecipeValue.c (revision 26747) +++ branches/eam_branches/psModules.stack.20100120/src/config/pmConfigRecipeValue.c (revision 26747) @@ -0,0 +1,84 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include + +#include "pmConfig.h" +#include "pmFPA.h" +#include "pmFPAview.h" + + +psMetadataItem *pmConfigRecipeValueByView(const pmConfig *config, // Configuration + const char *recipeName, // Name of recipe + const char *valueName, // Name of value in recipe + const pmFPA *fpa, // FPA of interest + const pmFPAview *view // View to component of interest + ) +{ + PS_ASSERT_PTR_NON_NULL(config, NULL); + PS_ASSERT_METADATA_NON_NULL(config->recipes, NULL); + PS_ASSERT_STRING_NON_EMPTY(recipeName, NULL); + PS_ASSERT_STRING_NON_EMPTY(valueName, NULL); + PS_ASSERT_PTR_NON_NULL(fpa, NULL); + PS_ASSERT_PTR_NON_NULL(view, NULL); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, recipeName); // Recipe of interest + if (!recipe) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", recipeName); + return NULL; + } + + psMetadataItem *fpaItem = psMetadataLookup(recipe, valueName); // Value for FPA or menu of chips + if (!fpaItem) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find value %s in recipe %s", valueName, recipeName); + return NULL; + } + if (view->chip == -1 && view->cell == -1) { + // Default value only + return fpaItem; + } + if (fpaItem->type != PS_DATA_METADATA) { + // Single value appropriate for all + return fpaItem; + } + + psMetadata *menu = fpaItem->data.md; // Menu with values by chip + + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip + + psMetadataItem *chipItem = psMetadataLookup(menu, chipName); // Value for chip of menu of cells + if (!chipItem) { + psError(PS_ERR_UNEXPECTED_NULL, true, + "Chip %s does not have a menu entry in %s within recipe %s.", + chipName, valueName, recipeName); + return NULL; + } + if (view->cell == -1) { + // Default value for chip + return chipItem; + } + if (chipItem->type != PS_DATA_METADATA) { + // Single value appropriate for all + return chipItem; + } + + menu = chipItem->data.md; // Menu with values by chip + + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + + psMetadataItem *cellItem = psMetadataLookup(menu, cellName); // Value for cell of menu of cells + if (!cellItem) { + psError(PS_ERR_UNEXPECTED_NULL, true, + "Chip %s, cell %s does not have a menu entry in %s within recipe %s --- using default.", + chipName, cellName, valueName, recipeName); + return NULL; + } + + return cellItem; +} Index: branches/eam_branches/psModules.stack.20100120/src/config/pmConfigRecipeValue.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/config/pmConfigRecipeValue.h (revision 26747) +++ branches/eam_branches/psModules.stack.20100120/src/config/pmConfigRecipeValue.h (revision 26747) @@ -0,0 +1,18 @@ +#ifndef PM_CONFIG_RECIPE_VALUE_H +#define PM_CONFIG_RECIPE_VALUE_H + +#include + +#include "pmConfig.h" +#include "pmFPA.h" +#include "pmFPAview.h" + +/// Return a recipe value according to the provided view (i.e., chip- and/or cell-dependent) +psMetadataItem *pmConfigRecipeValueByView(const pmConfig *config, // Configuration + const char *recipeName, // Name of recipe + const char *valueName, // Name of value in recipe + const pmFPA *fpa, // FPA of interest + const pmFPAview *view // View to component of interest + ); + +#endif Index: branches/eam_branches/psModules.stack.20100120/src/detrend/pmPattern.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/detrend/pmPattern.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/detrend/pmPattern.c (revision 26747) @@ -63,4 +63,5 @@ float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data + float background = stats->robustMedian; psFree(stats); psFree(rng); @@ -113,5 +114,5 @@ for (int x = 0; x < numCols; x++) { - image->data.F32[y][x] -= solution->data.F32[x]; + image->data.F32[y][x] += (background - solution->data.F32[x]); } psFree(solution); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmStackReject.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmStackReject.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmStackReject.c (revision 26747) @@ -35,6 +35,6 @@ { int size = kernels->size; // Half-size of convolution kernel - psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows, - xMin + size + 1, yMin + size + 1); // Polynomial + psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + size + 1, + yMin + size + 1); // Polynomial int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues, false, poorFrac); // Radius of bad box psTrace("psModules.imcombine", 10, "Growing by %d", box); @@ -150,4 +150,5 @@ pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image inRO->image = image; + convRO->image = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); for (int i = 0; i < numRegions; i++) { psRegion *region = subRegions->data[i]; // Region of interest @@ -165,8 +166,8 @@ // Image of the kernel at the centre of the region - float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) / - (float)kernels->numCols; - float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) / - (float)kernels->numRows; + float xNorm, yNorm; // Normalised coordinates + p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, 0.5 * (region->x1 - region->x0), + 0.5 * (region->y1 - region->y0), + kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax); psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false); if (!kernel) { Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.c (revision 26747) @@ -424,4 +424,6 @@ } else if (*source & subConvPoor) { *target |= maskPoor; + } else { + *target &= ~maskBad & ~maskPoor; } } @@ -574,19 +576,21 @@ } +void p_pmSubtractionPolynomialNormCoords(float *xOut, float *yOut, float xIn, float yIn, + int xMin, int xMax, int yMin, int yMax) +{ + float xNormSize = xMax - xMin, yNormSize = yMax - yMin; // Size to use for normalisation + *xOut = 2.0 * (float)(xIn - xMin - xNormSize/2.0) / xNormSize; + *yOut = 2.0 * (float)(yIn - yMin - yNormSize/2.0) / yNormSize; + return; +} + psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels, - int numCols, int numRows, int x, int y) + int x, int y) { assert(kernels); - assert(numCols > 0 && numRows > 0); - - // Size to use when calculating normalised coordinates (different from actual size when convolving - // subimage) - int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols); - int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows); - - // Normalised coordinates - float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize; - float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize; - + + float xNorm, yNorm; // Normalised coordinates + p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, x, y, + kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax); return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm); } @@ -1072,7 +1076,6 @@ // Only generate polynomial values every kernel footprint, since we have already assumed // (with the stamps) that it does not vary rapidly on this scale. - psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows, - xMin + x0 + size + 1, - yMin + y0 + size + 1); + psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + x0 + size + 1, + yMin + y0 + size + 1); // Polynomial float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term @@ -1158,4 +1161,5 @@ PM_ASSERT_READOUT_NON_NULL(ro1, false); PM_ASSERT_READOUT_IMAGE(ro1, false); + PM_ASSERT_READOUT_IMAGE(out1, false); numCols = ro1->image->numCols; numRows = ro1->image->numRows; @@ -1167,4 +1171,5 @@ PM_ASSERT_READOUT_NON_NULL(ro2, false); PM_ASSERT_READOUT_IMAGE(ro2, false); + PM_ASSERT_READOUT_IMAGE(out2, false); if (numCols == 0 && numRows == 0) { numCols = ro2->image->numCols; @@ -1200,41 +1205,10 @@ bool threaded = pmSubtractionThreaded(); // Running threaded? - // Outputs - if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out1->image) { - out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - if (ro1->variance) { - if (!out1->variance) { - out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - psImageInit(out1->variance, 0.0); - } - } - if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out2->image) { - out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - if (ro2->variance) { - if (!out2->variance) { - out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - psImageInit(out2->variance, 0.0); - } - } psImage *convMask = NULL; // Convolved mask image (common to inputs 1 and 2) if (subMask) { if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out1->mask) { - out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - } - psImageInit(out1->mask, 0); convMask = out1->mask; } if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - if (!out2->mask) { - out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - } - psImageInit(out2->mask, 0); if (!convMask) { convMask = out2->mask; @@ -1256,6 +1230,6 @@ // Get region for convolution: [xMin:xMax,yMin:yMax] - int xMin = size, xMax = numCols - size; - int yMin = size, yMax = numRows - size; + int xMin = kernels->xMin + size, xMax = kernels->xMax - size; + int yMin = kernels->yMin + size, yMax = kernels->yMax - size; if (region) { xMin = PS_MAX(region->x0, xMin); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.h (revision 26747) @@ -127,4 +127,11 @@ ); +/// Return normalised coordinates +void p_pmSubtractionPolynomialNormCoords( + float *xOut, float *yOut, ///< Normalised coordinates, returned + float xIn, float yIn, ///< Input coordinates + int xMin, int xMax, int yMin, int yMax ///< Bounds of validity + ); + /// Given (normalised) coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j psImage *p_pmSubtractionPolynomial(psImage *output, ///< Output matrix, or NULL @@ -138,5 +145,4 @@ psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, ///< Output matrix, or NULL const pmSubtractionKernels *kernels, ///< Kernel parameters - int numCols, int numRows, ///< Size of image of interest int x, int y ///< Position of interest ); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionAnalysis.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionAnalysis.c (revision 26747) @@ -118,5 +118,5 @@ // sample difference images - { + { psMetadataAddArray(analysis, PS_LIST_TAIL, "SUBTRACTION.SAMPLE.STAMP.SET", PS_META_DUPLICATE_OK, "Sample Difference Stamps", kernels->sampleStamps); } @@ -273,6 +273,5 @@ // Difference in background { - psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows, - numCols / 2.0, numRows / 2.0); // Polynomial + psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0); // Polynomial float bg = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background difference Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.c (revision 26747) @@ -28,4 +28,5 @@ static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated psVector *vector, // Least-squares vector, updated + double *norm, // Normalisation, updated const psKernel *input, // Input image (target) const psKernel *reference, // Reference image (convolution source) @@ -36,4 +37,5 @@ const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp + int normWindow, // Window (half-)size for normalisation measurement const pmSubtractionEquationCalculationMode mode ) @@ -156,7 +158,5 @@ if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { vector->data.F64[iIndex] = sumIC * poly[iTerm]; - // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B), - // instead, calculate A - \sum(k x B), with full hermitians - if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) { + if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { // subtract norm * sumRC * poly[iTerm] psAssert (kernels->solution1, "programming error: define solution first!"); @@ -174,4 +174,5 @@ double sumR = 0.0; // Sum of the reference double sumI = 0.0; // Sum of the input + double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -181,4 +182,10 @@ double rr = PS_SQR(ref); double one = 1.0; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { + normI1 += ref; + normI2 += in; + } + if (weight) { float wtVal = weight->kernel[y][x]; @@ -204,4 +211,7 @@ } } + + *norm = normI2 / normI1; + if (mode & PM_SUBTRACTION_EQUATION_NORM) { matrix->data.F64[normIndex][normIndex] = sumRR; @@ -217,4 +227,21 @@ matrix->data.F64[bgIndex][normIndex] = sumR; } + + // check for any NAN values in the result, skip if found: + for (int iy = 0; iy < matrix->numRows; iy++) { + for (int ix = 0; ix < matrix->numCols; ix++) { + if (!isfinite(matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + return false; + } + } + } + for (int ix = 0; ix < vector->n; ix++) { + if (!isfinite(vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + return false; + } + } + return true; } @@ -224,4 +251,5 @@ static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated psVector *vector, // Least-squares vector, updated + double *norm, // Normalisation, updated const psKernel *image1, // Image 1 const psKernel *image2, // Image 2 @@ -232,5 +260,7 @@ const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint // (Half-)Size of stamp + int footprint, // (Half-)Size of stamp + int normWindow, // Window (half-)size for normalisation measurement + const pmSubtractionEquationCalculationMode mode ) { @@ -272,4 +302,12 @@ } + + // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we + // choose to calculate + psImageInit(matrix, 0.0); + psVectorInit(vector, 1.0); + for (int i = 0; i < matrix->numCols; i++) { + matrix->data.F64[i][i] = 1.0; + } for (int i = 0; i < numKernels; i++) { @@ -306,19 +344,21 @@ } - // Spatial variation - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double aa = sumAA * poly2[iTerm][jTerm]; - double bb = sumBB * poly2[iTerm][jTerm]; - double ab = sumAB * poly2[iTerm][jTerm]; - - matrix->data.F64[iIndex][jIndex] = aa; - matrix->data.F64[jIndex][iIndex] = aa; - - matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; - matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; - - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; + // Spatial variation of kernel coeffs + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double aa = sumAA * poly2[iTerm][jTerm]; + double bb = sumBB * poly2[iTerm][jTerm]; + double ab = sumAB * poly2[iTerm][jTerm]; + + matrix->data.F64[iIndex][jIndex] = aa; + matrix->data.F64[jIndex][iIndex] = aa; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; + matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; + + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } } } @@ -340,10 +380,12 @@ } - // Spatial variation - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double ab = sumAB * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; + // Spatial variation of kernel coeffs + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double ab = sumAB * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } } } @@ -403,18 +445,32 @@ double bi2 = sumBI2 * poly[iTerm]; double ai1 = sumAI1 * poly[iTerm]; - double a = sumA * poly[iTerm]; + double a = sumA * poly[iTerm]; double bi1 = sumBI1 * poly[iTerm]; - double b = sumB * poly[iTerm]; - - matrix->data.F64[iIndex][normIndex] = ai1; - matrix->data.F64[normIndex][iIndex] = ai1; - matrix->data.F64[iIndex][bgIndex] = a; - matrix->data.F64[bgIndex][iIndex] = a; - matrix->data.F64[iIndex + numParams][normIndex] = bi1; - matrix->data.F64[normIndex][iIndex + numParams] = bi1; - matrix->data.F64[iIndex + numParams][bgIndex] = b; - matrix->data.F64[bgIndex][iIndex + numParams] = b; - vector->data.F64[iIndex] = ai2; - vector->data.F64[iIndex + numParams] = bi2; + double b = sumB * poly[iTerm]; + + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][normIndex] = ai1; + matrix->data.F64[normIndex][iIndex] = ai1; + matrix->data.F64[iIndex + numParams][normIndex] = bi1; + matrix->data.F64[normIndex][iIndex + numParams] = bi1; + } + if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { + matrix->data.F64[iIndex][bgIndex] = a; + matrix->data.F64[bgIndex][iIndex] = a; + matrix->data.F64[iIndex + numParams][bgIndex] = b; + matrix->data.F64[bgIndex][iIndex + numParams] = b; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + vector->data.F64[iIndex] = ai2; + vector->data.F64[iIndex + numParams] = bi2; + if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { + // subtract norm * sumRC * poly[iTerm] + psAssert (kernels->solution1, "programming error: define solution first!"); + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation + vector->data.F64[iIndex] -= norm * ai1; + vector->data.F64[iIndex + numParams] -= norm * bi1; + } + } } } @@ -425,4 +481,5 @@ double sumI2 = 0.0; // Sum of I_2 (for vector, background) double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) + double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -433,4 +490,9 @@ double one = 1.0; double i1i2 = i1 * i2; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { + normI1 += i1; + normI2 += i2; + } if (weight) { @@ -457,10 +519,36 @@ } } - matrix->data.F64[bgIndex][normIndex] = sumI1; - matrix->data.F64[normIndex][bgIndex] = sumI1; - matrix->data.F64[normIndex][normIndex] = sumI1I1; - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI2; - vector->data.F64[normIndex] = sumI1I2; + + *norm = normI2 / normI1; + + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + matrix->data.F64[normIndex][normIndex] = sumI1I1; + vector->data.F64[normIndex] = sumI1I2; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + matrix->data.F64[bgIndex][bgIndex] = sum1; + vector->data.F64[bgIndex] = sumI2; + } + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { + matrix->data.F64[bgIndex][normIndex] = sumI1; + matrix->data.F64[normIndex][bgIndex] = sumI1; + } + + // check for any NAN values in the result, skip if found: + for (int iy = 0; iy < matrix->numRows; iy++) { + for (int ix = 0; ix < matrix->numCols; ix++) { + if (!isfinite(matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + return false; + } + } + } + for (int ix = 0; ix < vector->n; ix++) { + if (!isfinite(vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + return false; + } + } + return true; @@ -469,5 +557,5 @@ #if 1 // Add in penalty term to least-squares vector -static bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term +bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term psVector *vector, // Vector to which to add in penalty term const pmSubtractionKernels *kernels, // Kernel parameters @@ -482,4 +570,17 @@ int spatialOrder = kernels->spatialOrder; // Order of spatial variations int numKernels = kernels->num; // Number of kernel components + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations + int numParams = numKernels * numSpatial; // Number of kernel parameters + + // order is : + // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0] + // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] + // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] + // [norm] + // [bg] + // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] + // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] + // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1] + for (int i = 0; i < numKernels; i++) { for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { @@ -488,4 +589,10 @@ psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty"); matrix->data.F64[index][index] += norm * penalties->data.F32[i]; + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i]; + // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) + // penalties scale with second moments + // + } } } @@ -668,17 +775,18 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: - status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1, + status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, weight, window, stamp->convolutions1, kernels, - polyValues, footprint, mode); + polyValues, footprint, stamps->normWindow, mode); break; case PM_SUBTRACTION_MODE_2: - status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2, + status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, weight, window, stamp->convolutions2, kernels, - polyValues, footprint, mode); + polyValues, footprint, stamps->normWindow, mode); break; case PM_SUBTRACTION_MODE_DUAL: - status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2, + status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, + stamp->image1, stamp->image2, weight, window, stamp->convolutions1, stamp->convolutions2, - kernels, polyValues, footprint); + kernels, polyValues, footprint, stamps->normWindow, mode); break; default: @@ -721,5 +829,6 @@ } -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode) +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, + const pmSubtractionEquationCalculationMode mode) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -843,349 +952,6 @@ psVectorInit(sumVector, 0.0); psImageInit(sumMatrix, 0.0); - int numStamps = 0; // Number of good stamps - for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - - if (stamp->status == PM_SUBTRACTION_STAMP_USED) { - -#ifdef TESTING - // XXX double-check for NAN in data: - for (int iy = 0; iy < stamp->matrix->numRows; iy++) { - for (int ix = 0; ix < stamp->matrix->numCols; ix++) { - if (!isfinite(stamp->matrix->data.F64[iy][ix])) { - fprintf (stderr, "WARNING: NAN in matrix\n"); - } - } - } - for (int ix = 0; ix < stamp->vector->n; ix++) { - if (!isfinite(stamp->vector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - - (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); - (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); - numStamps++; - } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); - } - } - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - -#if 0 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); -#endif - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } - { - psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); - psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); - psFree(inverse); - } - { - psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); - psImage *Xt = psMatrixTranspose(NULL, X); - psImage *XtX = psMatrixMultiply(NULL, Xt, X); - psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); - psFree(X); - psFree(Xt); - psFree(XtX); - } -#endif - -# define SVD_ANALYSIS 0 -# define COEFF_SIG 0.0 -# define SVD_TOL 0.0 - - psVector *solution = NULL; - psVector *solutionErr = NULL; - -#ifdef TESTING - psFitsWriteImageSimple("A.fits", sumMatrix, NULL); - psVectorWriteFile ("B.dat", sumVector); -#endif - - // Use SVD to determine the kernel coeffs (and validate) - if (SVD_ANALYSIS) { - - // We have sumVector and sumMatrix. we are trying to solve the following equation: - // sumMatrix * x = sumVector. - - // we can use any standard matrix inversion to solve this. However, the basis - // functions in general have substantial correlation, so that the solution may be - // somewhat poorly determined or unstable. If not numerically ill-conditioned, the - // system of equations may be statistically ill-conditioned. Noise in the image - // will drive insignificant, but correlated, terms in the solution. To avoid these - // problems, we can use SVD to identify numerically unconstrained values and to - // avoid statistically badly determined value. - - // A = sumMatrix, B = sumVector - // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T - // x = V (1/w) (U^T B) - // \sigma_x = sqrt(diag(A^{-1})) - // solve for x and A^{-1} to get x & dx - // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 - // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 - - // If I use the SVD trick to re-condition the matrix, I need to break out the - // kernel and normalization terms from the background term. - // XXX is this true? or was this due to an error in the analysis? - - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - - // now pull out the kernel elements into their own square matrix - psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); - psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { - if (ix == bgIndex) continue; - for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { - if (iy == bgIndex) continue; - kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; - ky++; - } - kernelVector->data.F64[kx] = sumVector->data.F64[ix]; - kx++; - } - - psImage *U = NULL; - psImage *V = NULL; - psVector *w = NULL; - if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { - psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); - return NULL; - } - - // calculate A_inverse: - // Ainv = V * w * U^T - psImage *wUt = p_pmSubSolve_wUt (w, U); - psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); - psImage *Xvar = NULL; - psFree (wUt); - -# ifdef TESTING - // kernel terms: - for (int i = 0; i < w->n; i++) { - fprintf (stderr, "w: %f\n", w->data.F64[i]); - } -# endif - // loop over w adding in more and more of the values until chisquare is no longer - // dropping significantly. - // XXX this does not seem to work very well: we seem to need all terms even for - // simple cases... - - psVector *Xsvd = NULL; - { - psVector *Ax = NULL; - psVector *UtB = NULL; - psVector *wUtB = NULL; - - psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); - psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); - psVectorInit (wMask, 1); // start by masking everything - - double chiSquareLast = NAN; - int maxWeight = 0; - - double Axx, Bx, y2; - - // XXX this is an attempt to exclude insignificant modes. - // it was not successful with the ISIS kernel set: removing even - // the least significant mode leaves additional ringing / noise - // because the terms are so coupled. - for (int k = 0; false && (k < w->n); k++) { - - // unmask the k-th weight - wMask->data.U8[k] = 0; - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; - chiSquareLast = chiSquare; - - // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); - if (k && !maxWeight && (deltaChi < 1.0)) { - maxWeight = k; - } - } - - // keep all terms or we get extra ringing - maxWeight = w->n; - psVectorInit (wMask, 1); - for (int k = 0; k < maxWeight; k++) { - wMask->data.U8[k] = 0; - } - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); - - // re-calculate A^{-1} to get new variances: - // Ainv = V * w * U^T - // XXX since we keep all terms, this is identical to Ainv - psImage *wUt = p_pmSubSolve_wUt (wApply, U); - Xvar = p_pmSubSolve_VwUt (V, wUt); - psFree (wUt); - - psFree (Ax); - psFree (UtB); - psFree (wUtB); - psFree (wApply); - psFree (wMask); - } - - // copy the kernel solutions to the full solution vector: - solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { - if (ix == bgIndex) { - solution->data.F64[ix] = 0; - solutionErr->data.F64[ix] = 0.001; - continue; - } - solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); - solution->data.F64[ix] = Xsvd->data.F64[kx]; - kx++; - } - - psFree (kernelMatrix); - psFree (kernelVector); - - psFree (U); - psFree (V); - psFree (w); - - psFree (Ainv); - psFree (Xsvd); - } else { - psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); - if (!luMatrix) { - psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - psFree(luMatrix); - psFree(permutation); - return NULL; - } - - solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); - psFree(luMatrix); - psFree(permutation); - if (!solution) { - psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - return NULL; - } - - // XXX LUD does not provide A^{-1}? fake the error for now - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - for (int ix = 0; ix < sumVector->n; ix++) { - solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; - } - } - - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); - psVectorInit (kernels->solution1, 0.0); - } - - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - int spatialOrder = kernels->spatialOrder; // Order of spatial variation - int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms - for (int i = 0; i < numKernels * numPoly; i++) { - // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); - if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { - // XXX fprintf (stderr, "drop\n"); - kernels->solution1->data.F64[i] = 0.0; - } else { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } - } - // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); - // fprintf (stderr, "chi square Brute = %f\n", chiSquare); - - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - -#ifdef TESTING - // XXX double-check for NAN in data: - for (int ix = 0; ix < kernels->solution1->n; ix++) { - if (!isfinite(kernels->solution1->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - - } else { - // Dual convolution solution - - // Accumulation of stamp matrices/vectors - psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); - psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); - psImageInit(sumMatrix, 0.0); - psVectorInit(sumVector, 0.0); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations int numStamps = 0; // Number of good stamps @@ -1195,4 +961,122 @@ (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + psVectorAppend(norms, stamp->norm); + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); + numStamps++; + } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); + } + } + +#if 0 + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); +#endif + + psVector *solution = NULL; // Solution to equation! + solution = psVectorAlloc(numParams, PS_TYPE_F64); + psVectorInit(solution, 0); + +#if 0 + // Regular, straight-forward solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); +#else + { + // Solve normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(sumMatrix); + psFree(sumVector); + psFree(norms); + return false; + } + + double normValue = stats->robustMedian; + // double bgValue = 0.0; + + psFree(stats); + + fprintf(stderr, "Norm: %lf\n", normValue); + + // Solve kernel components + for (int i = 0; i < numSolution1; i++) { + sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; + + sumMatrix->data.F64[i][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i] = 0.0; + } + sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; + sumMatrix->data.F64[bgIndex][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][bgIndex] = 0.0; + + sumMatrix->data.F64[normIndex][normIndex] = 1.0; + sumVector->data.F64[normIndex] = 0.0; + + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + + solution->data.F64[normIndex] = normValue; + } +# endif + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64); + psVectorInit(kernels->solution1, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int ix = 0; ix < kernels->solution1->n; ix++) { + if (!isfinite(kernels->solution1->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + + } else { + // Dual convolution solution + + // Accumulation of stamp matrices/vectors + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix, 0.0); + psVectorInit(sumVector, 0.0); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + + int numStamps = 0; // Number of good stamps + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + + psVectorAppend(norms, stamp->norm); pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); @@ -1207,6 +1091,9 @@ #if 1 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); + // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); + + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0); #endif @@ -1216,74 +1103,74 @@ #if 0 + // Regular, straight-forward solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); +#else { - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - if (!solution) { - psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n"); + // Solve normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + +#if 0 + psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); + psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); + + normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; + normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; + normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; + + normVector->data.F64[0] = sumVector->data.F64[normIndex]; + normVector->data.F64[1] = sumVector->data.F64[bgIndex]; + + psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); + + double normValue = normSolution->data.F64[0]; + double bgValue = normSolution->data.F64[1]; + + psFree(normMatrix); + psFree(normVector); + psFree(normSolution); +#endif + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation"); + psFree(stats); psFree(sumMatrix); psFree(sumVector); - return NULL; - } - -#ifdef TESTING - for (int i = 0; i < numParams; i++) { - fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]); - } -#endif - - int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations - int numKernels = kernels->num; // Number of kernel basis functions - -#define MASK_BASIS(INDEX) \ - { \ - for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \ - for (int k = 0; k < numParams; k++) { \ - sumMatrix->data.F64[index][k] = sumMatrix->data.F64[k][index] = 0.0; \ - } \ - sumMatrix->data.F64[index][index] = 1.0; \ - sumVector->data.F64[index] = 0.0; \ - } \ - } - - #define TOL 1.0e-3 - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(solution->data.F64[normIndex]); // Normalisation - double thresh = norm * TOL; // Threshold for low parameters - for (int i = 0; i < numKernels; i++) { - // Getting 0th order parameter value. In the presence of spatial variation, the actual value - // of the parameter will vary over the image. We are in effect getting the value in the - // centre of the image. If we use different polynomial functions (e.g., Chebyshev), we may - // have to change this to properly determine the value of the parameter at the centre. - double param1 = solution->data.F64[i], - param2 = solution->data.F64[numSolution1 + i]; // Parameters of interest - bool mask1 = false, mask2 = false; // Masked the parameter? - if (fabs(param1) < thresh) { - psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i); - MASK_BASIS(i); - mask1 = true; - } - if (fabs(param2) < thresh) { - psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i); - MASK_BASIS(numSolution1 + i); - mask2 = true; - } - - if (!mask1 && !mask2) { - if (fabs(param1) > fabs(param2)) { - psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i); - MASK_BASIS(numSolution1 + i); - } else { - psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i); - MASK_BASIS(i); - } - } - } - -#ifdef TESTING - psFitsWriteImageSimple ("sumMatrixFix.fits", sumMatrix, NULL); - psVectorWriteFile("sumVectorFix.dat", sumVector); -#endif - } -#endif /*** kernel-masking block ***/ - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + psFree(norms); + return false; + } + + double normValue = stats->robustMedian; + + psFree(stats); + + fprintf(stderr, "Norm: %lf\n", normValue); + + // Solve kernel components + for (int i = 0; i < numSolution2; i++) { + sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; + sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; + + sumMatrix->data.F64[i][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i] = 0.0; + + sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; + } + sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; + sumMatrix->data.F64[bgIndex][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][bgIndex] = 0.0; + + sumMatrix->data.F64[normIndex][normIndex] = 1.0; + + sumVector->data.F64[normIndex] = 0.0; + + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + + solution->data.F64[normIndex] = normValue; + } +#endif + #ifdef TESTING @@ -1296,10 +1183,33 @@ psFree(sumVector); + psFree(norms); + if (!kernels->solution1) { kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64); + psVectorInit (kernels->solution1, 0.0); } if (!kernels->solution2) { kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64); - } + psVectorInit (kernels->solution2, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + for (int i = 0; i < numKernels * numSpatial; i++) { + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; + } + } + memcpy(kernels->solution1->data.F64, solution->data.F64, @@ -1364,4 +1274,9 @@ } float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix)); + if (!isfinite(sum)) return false; + if (!isfinite(dmax)) return false; + if (!isfinite(dmin)) return false; + if (!isfinite(peak)) return false; + // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak); psVectorAppend(fSigRes, sigma/sum); @@ -1973,2 +1888,311 @@ } + +# if 0 + +#ifdef TESTING + psFitsWriteImageSimple("A.fits", sumMatrix, NULL); + psVectorWriteFile ("B.dat", sumVector); +#endif + +# define SVD_ANALYSIS 0 +# define COEFF_SIG 0.0 +# define SVD_TOL 0.0 + + // Use SVD to determine the kernel coeffs (and validate) + if (SVD_ANALYSIS) { + + // We have sumVector and sumMatrix. we are trying to solve the following equation: + // sumMatrix * x = sumVector. + + // we can use any standard matrix inversion to solve this. However, the basis + // functions in general have substantial correlation, so that the solution may be + // somewhat poorly determined or unstable. If not numerically ill-conditioned, the + // system of equations may be statistically ill-conditioned. Noise in the image + // will drive insignificant, but correlated, terms in the solution. To avoid these + // problems, we can use SVD to identify numerically unconstrained values and to + // avoid statistically badly determined value. + + // A = sumMatrix, B = sumVector + // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T + // x = V (1/w) (U^T B) + // \sigma_x = sqrt(diag(A^{-1})) + // solve for x and A^{-1} to get x & dx + // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 + // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 + + // If I use the SVD trick to re-condition the matrix, I need to break out the + // kernel and normalization terms from the background term. + // XXX is this true? or was this due to an error in the analysis? + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + + // now pull out the kernel elements into their own square matrix + psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); + psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { + if (ix == bgIndex) continue; + for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { + if (iy == bgIndex) continue; + kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; + ky++; + } + kernelVector->data.F64[kx] = sumVector->data.F64[ix]; + kx++; + } + + psImage *U = NULL; + psImage *V = NULL; + psVector *w = NULL; + if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { + psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); + return NULL; + } + + // calculate A_inverse: + // Ainv = V * w * U^T + psImage *wUt = p_pmSubSolve_wUt (w, U); + psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); + psImage *Xvar = NULL; + psFree (wUt); + +# ifdef TESTING + // kernel terms: + for (int i = 0; i < w->n; i++) { + fprintf (stderr, "w: %f\n", w->data.F64[i]); + } +# endif + // loop over w adding in more and more of the values until chisquare is no longer + // dropping significantly. + // XXX this does not seem to work very well: we seem to need all terms even for + // simple cases... + + psVector *Xsvd = NULL; + { + psVector *Ax = NULL; + psVector *UtB = NULL; + psVector *wUtB = NULL; + + psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); + psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); + psVectorInit (wMask, 1); // start by masking everything + + double chiSquareLast = NAN; + int maxWeight = 0; + + double Axx, Bx, y2; + + // XXX this is an attempt to exclude insignificant modes. + // it was not successful with the ISIS kernel set: removing even + // the least significant mode leaves additional ringing / noise + // because the terms are so coupled. + for (int k = 0; false && (k < w->n); k++) { + + // unmask the k-th weight + wMask->data.U8[k] = 0; + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; + chiSquareLast = chiSquare; + + // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); + if (k && !maxWeight && (deltaChi < 1.0)) { + maxWeight = k; + } + } + + // keep all terms or we get extra ringing + maxWeight = w->n; + psVectorInit (wMask, 1); + for (int k = 0; k < maxWeight; k++) { + wMask->data.U8[k] = 0; + } + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); + + // re-calculate A^{-1} to get new variances: + // Ainv = V * w * U^T + // XXX since we keep all terms, this is identical to Ainv + psImage *wUt = p_pmSubSolve_wUt (wApply, U); + Xvar = p_pmSubSolve_VwUt (V, wUt); + psFree (wUt); + + psFree (Ax); + psFree (UtB); + psFree (wUtB); + psFree (wApply); + psFree (wMask); + } + + // copy the kernel solutions to the full solution vector: + solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { + if (ix == bgIndex) { + solution->data.F64[ix] = 0; + solutionErr->data.F64[ix] = 0.001; + continue; + } + solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); + solution->data.F64[ix] = Xsvd->data.F64[kx]; + kx++; + } + + psFree (kernelMatrix); + psFree (kernelVector); + + psFree (U); + psFree (V); + psFree (w); + + psFree (Ainv); + psFree (Xsvd); + } else { + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); + if (!luMatrix) { + psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + + solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); + psFree(luMatrix); + psFree(permutation); + if (!solution) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + return NULL; + } + + // XXX LUD does not provide A^{-1}? fake the error for now + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + for (int ix = 0; ix < sumVector->n; ix++) { + solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; + } + } + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); + psVectorInit (kernels->solution1, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); + if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { + // XXX fprintf (stderr, "drop\n"); + kernels->solution1->data.F64[i] = 0.0; + } else { + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + } + // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); + // fprintf (stderr, "chi square Brute = %f\n", chiSquare); + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); +# endif + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int iy = 0; iy < stamp->matrix->numRows; iy++) { + for (int ix = 0; ix < stamp->matrix->numCols; ix++) { + if (!isfinite(stamp->matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + } + } + } + for (int ix = 0; ix < stamp->vector->n; ix++) { + if (!isfinite(stamp->vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + +#ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + +#ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } + { + psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); + psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); + psFree(inverse); + } + { + psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); + psImage *Xt = psMatrixTranspose(NULL, X); + psImage *XtX = psMatrixMultiply(NULL, Xt, X); + psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); + psFree(X); + psFree(Xt); + psFree(XtX); + } +#endif + Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.h (revision 26747) @@ -10,5 +10,5 @@ PM_SUBTRACTION_EQUATION_BG = 0x02, PM_SUBTRACTION_EQUATION_KERNELS = 0x04, - PM_SUBTRACTION_EQUATION_ALL = 0x05, // value should be NORM | BG | KERNELS + PM_SUBTRACTION_EQUATION_ALL = 0x07, // value should be NORM | BG | KERNELS } pmSubtractionEquationCalculationMode; @@ -21,5 +21,5 @@ const pmSubtractionKernels *kernels, ///< Kernel parameters int index, ///< Index of stamp - const pmSubtractionEquationCalculationMode mode + const pmSubtractionEquationCalculationMode mode ); @@ -27,5 +27,5 @@ bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps const pmSubtractionKernels *kernels, ///< Kernel parameters - const pmSubtractionEquationCalculationMode mode + const pmSubtractionEquationCalculationMode mode ); @@ -33,5 +33,5 @@ bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters const pmSubtractionStampList *stamps, ///< Stamps - const pmSubtractionEquationCalculationMode mode + const pmSubtractionEquationCalculationMode mode ); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionIO.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionIO.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionIO.c (revision 26747) @@ -80,10 +80,5 @@ 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); } @@ -126,6 +121,4 @@ 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_2) { psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL1, 0, "Solution vector 1", kernel->solution1); @@ -318,6 +311,4 @@ TABLE_LOOKUP(int, S32, bg, NAME_BG); - TABLE_LOOKUP(int, S32, numCols, NAME_COLS); - TABLE_LOOKUP(int, S32, numRows, NAME_ROWS); TABLE_LOOKUP(float, F32, mean, NAME_MEAN); @@ -325,7 +316,5 @@ TABLE_LOOKUP(int, S32, numStamps, NAME_NUMSTAMPS); - pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode); - kernels->numCols = numCols; - kernels->numRows = numRows; + pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, *region, mode); kernels->mean = mean; kernels->rms = rms; Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.c (revision 26747) @@ -276,5 +276,5 @@ pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, const psVector *fwhmsIN, const psVector *ordersIN, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); @@ -306,5 +306,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); @@ -332,5 +332,5 @@ pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder, const psVector *fwhmsIN, const psVector *ordersIN, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); @@ -362,5 +362,6 @@ } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, spatialOrder, penalty, mode); // The kernels + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, + spatialOrder, penalty, bounds, mode); // Kernels psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); @@ -389,5 +390,5 @@ pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder, const psVector *fwhmsIN, const psVector *ordersIN, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); @@ -418,8 +419,10 @@ } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, spatialOrder, penalty, mode); // The kernels + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, + spatialOrder, penalty, bounds, mode); // Kernels psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); - psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", params, spatialOrder, num); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", + params, spatialOrder, num); psFree(params); @@ -440,6 +443,6 @@ pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder, - const psVector *fwhmsIN, const psVector *ordersIN, - float penalty, pmSubtractionMode mode) + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); @@ -470,5 +473,6 @@ } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, spatialOrder, penalty, mode); // The kernels + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, + spatialOrder, penalty, bounds, mode); // Kernels psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); @@ -545,5 +549,5 @@ pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, - int size, int spatialOrder, float penalty, + int size, int spatialOrder, float penalty, psRegion bounds, pmSubtractionMode mode) { @@ -559,4 +563,8 @@ kernels->uStop = NULL; kernels->vStop = NULL; + kernels->xMin = bounds.x0; + kernels->xMax = bounds.x1; + kernels->yMin = bounds.y0; + kernels->yMax = bounds.y1; kernels->preCalc = psArrayAlloc(numBasisFunctions); kernels->penalty = penalty; @@ -567,6 +575,4 @@ kernels->bgOrder = 0; kernels->mode = mode; - kernels->numCols = 0; - kernels->numRows = 0; kernels->solution1 = NULL; kernels->solution2 = NULL; @@ -641,5 +647,5 @@ } -pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, psRegion bounds, pmSubtractionMode mode) { @@ -650,5 +656,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty); psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", @@ -665,8 +671,8 @@ pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, const psVector *fwhms, const psVector *orders, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, - penalty, mode); // Kernels + penalty, bounds, mode); // Kernels if (!kernels) { return NULL; @@ -677,5 +683,5 @@ pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -698,5 +704,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder, @@ -769,5 +775,5 @@ pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty, - pmSubtractionMode mode) + psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -796,5 +802,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty); @@ -865,5 +871,5 @@ pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, const psVector *orders, int inner, float penalty, - pmSubtractionMode mode) + psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -878,5 +884,5 @@ pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, - penalty, mode); // Kernels + penalty, bounds, mode); // Kernels kernels->type = PM_SUBTRACTION_KERNEL_GUNK; psStringPrepend(&kernels->description, "GUNK="); @@ -894,5 +900,5 @@ // RINGS --- just what it says pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -925,5 +931,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder, @@ -1047,26 +1053,26 @@ pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *fwhms, const psVector *orders, int inner, - int binning, int ringsOrder, float penalty, + int binning, int ringsOrder, float penalty, psRegion bounds, pmSubtractionMode mode) { switch (type) { case PM_SUBTRACTION_KERNEL_POIS: - return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode); + return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_ISIS: - return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode); + return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: - return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, mode); + return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_HERM: - return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode); + return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_DECONV_HERM: - return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode); + return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_SPAM: - return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_FRIES: - return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode); + return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_GUNK: - return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode); + return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_RINGS: - return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode); + return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); @@ -1104,5 +1110,5 @@ pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder, - pmSubtractionMode mode) + psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_STRING_NON_EMPTY(description, NULL); @@ -1129,5 +1135,5 @@ type = pmSubtractionKernelsTypeFromString (description); - char *ptr = strchr(description, '('); + char *ptr = strchr(description, '(') + 1; psAssert (ptr, "description is missing kernel parameters"); @@ -1157,7 +1163,5 @@ PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); penalty = parseStringFloat(ptr); - - return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); - + break; case PM_SUBTRACTION_KERNEL_RINGS: PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); @@ -1166,9 +1170,11 @@ PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); - return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); + break; default: psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported."); } - return NULL; + + return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, + ringsOrder, penalty, bounds, mode); } @@ -1246,6 +1252,8 @@ out->bgOrder = in->bgOrder; out->mode = in->mode; - out->numCols = in->numCols; - out->numRows = in->numRows; + out->xMin = in->xMin; + out->xMax = in->xMax; + out->yMin = in->yMin; + out->yMax = in->yMax; out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL; out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL; Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.h (revision 26747) @@ -12,5 +12,5 @@ PM_SUBTRACTION_KERNEL_ISIS_RADIAL, ///< ISIS + higher-order radial Hermitians PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels - PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved Hermitian polynomial kernels + PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved Hermitian polynomial kernels PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction @@ -32,4 +32,5 @@ pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels psString description; ///< Description of the kernel parameters + int xMin, xMax, yMin, yMax; ///< Bounds of image (for normalisation) long num; ///< Number of kernel components (not including the spatial ones) psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM) @@ -44,27 +45,26 @@ int bgOrder; ///< The order for the background fitting pmSubtractionMode mode; ///< Mode for subtraction - int numCols, numRows; ///< Size of image (for normalisation), or zero to use image provided psVector *solution1, *solution2; ///< Solution for the PSF matching // Quality information float mean, rms; ///< Mean and RMS of chi^2 from stamps int numStamps; ///< Number of good stamps - float fSigResMean; ///< mean fractional stdev of residuals - float fSigResStdev; ///< stdev of fractional stdev of residuals - float fMaxResMean; ///< mean fractional positive swing in residuals - float fMaxResStdev; ///< stdev of fractional positive swing in residuals - float fMinResMean; ///< mean fractional negative swing in residuals - float fMinResStdev; ///< stdev of fractional negative swing in residuals - psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations + float fSigResMean; ///< mean fractional stdev of residuals + float fSigResStdev; ///< stdev of fractional stdev of residuals + float fMaxResMean; ///< mean fractional positive swing in residuals + float fMaxResStdev; ///< stdev of fractional positive swing in residuals + float fMinResMean; ///< mean fractional negative swing in residuals + float fMinResStdev; ///< stdev of fractional negative swing in residuals + psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations } pmSubtractionKernels; // pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures typedef struct { - psVector *uCoords; // used by RINGS - psVector *vCoords; // used by RINGS - psVector *poly; // used by RINGS - - psVector *xKernel; // used by ISIS, HERM, DECONV_HERM - psVector *yKernel; // used by ISIS, HERM, DECONV_HERM - psKernel *kernel; // used by ISIS, HERM, DECONV_HERM + psVector *uCoords; // used by RINGS + psVector *vCoords; // used by RINGS + psVector *poly; // used by RINGS + + psVector *xKernel; // used by ISIS, HERM, DECONV_HERM + psVector *yKernel; // used by ISIS, HERM, DECONV_HERM + psKernel *kernel; // used by ISIS, HERM, DECONV_HERM } pmSubtractionKernelPreCalc; @@ -162,4 +162,5 @@ int spatialOrder, ///< Order of spatial variations float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -168,8 +169,8 @@ pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc( pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated) - int uOrder, ///< order in x-direction - int vOrder, ///< order in x-direction - int size, ///< Half-size of the kernel - float sigma ///< sigma of gaussian kernel + int uOrder, ///< order in x-direction + int vOrder, ///< order in x-direction + int size, ///< Half-size of the kernel + float sigma ///< sigma of gaussian kernel ); @@ -179,4 +180,5 @@ int spatialOrder, ///< Order of spatial variations float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -188,4 +190,5 @@ const psVector *orders, ///< Polynomial order of gaussians float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -197,4 +200,5 @@ const psVector *orders, ///< Polynomial order of gaussians float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -202,9 +206,10 @@ /// Generate ISIS + RADIAL_HERM kernels pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - const psVector *fwhms, ///< Gaussian FWHMs - const psVector *orders, ///< Polynomial order of gaussians - float penalty, ///< Penalty for wideness - pmSubtractionMode mode ///< Mode for subtraction + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< Polynomial order of gaussians + float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity + pmSubtractionMode mode ///< Mode for subtraction ); @@ -215,4 +220,5 @@ const psVector *orders, ///< order of hermitian polynomials float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -220,9 +226,10 @@ /// Generate DECONV_HERM kernels pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - const psVector *fwhms, ///< Gaussian FWHMs - const psVector *orders, ///< order of hermitian polynomials - float penalty, ///< Penalty for wideness - pmSubtractionMode mode ///< Mode for subtraction + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< order of hermitian polynomials + float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity + pmSubtractionMode mode ///< Mode for subtraction ); @@ -233,4 +240,5 @@ int binning, ///< Kernel binning factor float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -241,4 +249,5 @@ int inner, ///< Inner radius to preserve unbinned float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -251,4 +260,5 @@ int inner, ///< Inner radius containing grid of delta functions float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -260,4 +270,5 @@ int ringsOrder, ///< Polynomial order float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -274,4 +285,5 @@ int ringsOrder, ///< Polynomial order for RINGS float penalty, ///< Penalty for wideness + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); @@ -281,4 +293,5 @@ const char *description, ///< Description of kernel int bgOrder, ///< Polynomial order for background fitting + psRegion bounds, ///< Bounds for validity pmSubtractionMode mode ///< Mode for subtraction ); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMask.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMask.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMask.c (revision 26747) @@ -38,14 +38,15 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psImageMaskType maskVal, - int size, int footprint, float badFrac, pmSubtractionMode mode) -{ - PS_ASSERT_IMAGE_NON_NULL(mask1, NULL); - PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_IMAGE_MASK, NULL); - if (mask2) { - PS_ASSERT_IMAGE_NON_NULL(mask2, NULL); - PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_IMAGE_MASK, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL); - } +psImage *pmSubtractionMask(psRegion *bounds, const pmReadout *ro1, const pmReadout *ro2, + psImageMaskType maskVal, int size, int footprint, float badFrac, + pmSubtractionMode mode) +{ + PM_ASSERT_READOUT_NON_NULL(ro1, NULL); + PM_ASSERT_READOUT_IMAGE(ro1, NULL); + PM_ASSERT_READOUT_MASK(ro1, NULL); + PM_ASSERT_READOUT_NON_NULL(ro2, NULL); + PM_ASSERT_READOUT_IMAGE(ro2, NULL); + PM_ASSERT_READOUT_MASK(ro2, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL); PS_ASSERT_INT_NONNEGATIVE(size, NULL); PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); @@ -55,28 +56,38 @@ } - int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of the images // Dereference inputs for convenience - psImageMaskType **data1 = mask1->data.PS_TYPE_IMAGE_MASK_DATA; - psImageMaskType **data2 = NULL; - if (mask2) { - data2 = mask2->data.PS_TYPE_IMAGE_MASK_DATA; - } + psF32 **imageData1 = ro1->image->data.F32, **imageData2 = ro2->image->data.F32; + psImageMaskType **maskData1 = ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA, + **maskData2 = ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA; // First, a pass through to determine the fraction of bad pixels - if (isfinite(badFrac) && badFrac != 1.0) { + if (bounds || (isfinite(badFrac) && badFrac != 1.0)) { + int xMin = numCols, xMax = 0, yMin = numRows, yMax = 0; // Bounds of good pixels int numBad = 0; // Number of bad pixels for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - if (data1[y][x] & maskVal) { + if ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x])) { numBad++; continue; } - if (data2 && data2[y][x] & maskVal) { + if ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x])) { numBad++; + continue; } - } - } - if (numBad > badFrac * numCols * numRows) { + xMin = PS_MIN(xMin, x); + xMax = PS_MAX(xMax, x); + yMin = PS_MIN(yMin, y); + yMax = PS_MAX(yMax, y); + } + } + if (bounds) { + bounds->x0 = xMin; + bounds->x1 = xMax; + bounds->y0 = yMin; + bounds->y1 = yMax; + } + if (isfinite(badFrac) && badFrac != 1.0 && numBad > badFrac * numCols * numRows) { psError(PM_ERR_SMALL_AREA, true, "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", @@ -117,8 +128,8 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - if (data1[y][x] & maskVal) { + if (maskData1[y][x] & maskVal) { maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1; } - if (data2 && data2[y][x] & maskVal) { + if (maskData2[y][x] & maskVal) { maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2; } Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMask.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMask.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMask.h (revision 26747) @@ -5,11 +5,13 @@ /// Generate a mask for use in the subtraction process -psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved) - const psImage *inMask, ///< Mask for the input image, or NULL - psImageMaskType maskVal, ///< Value to mask out - int size, ///< Half-size of the kernel (pmSubtractionKernels.size) - int footprint, ///< Half-size of the kernel footprint - float badFrac, ///< Maximum fraction of bad input pixels to accept - pmSubtractionMode mode ///< Subtraction mode +psImage *pmSubtractionMask( + psRegion *bounds, ///< Bounds of valid pixels (or NULL), returned + const pmReadout *ro1, ///< Readout 1 + const pmReadout *ro2, ///< Readout 2 + psImageMaskType maskVal, ///< Value to mask out + int size, ///< Half-size of the kernel (pmSubtractionKernels.size) + int footprint, ///< Half-size of the kernel footprint + float badFrac, ///< Maximum fraction of bad input pixels to accept + pmSubtractionMode mode ///< Subtraction mode ); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.c (revision 26747) @@ -33,5 +33,4 @@ # else # define SUBMODE PM_SUBTRACTION_EQUATION_ALL -// # define SUBMODE PM_SUBTRACTION_EQUATION_KERNELS # endif @@ -71,8 +70,9 @@ const psImage *subMask, // Mask for subtraction, or NULL psImage *variance, // Variance map - const psRegion *region, // Region of interest, or NULL + const psRegion *region, // Region of interest float thresh1, // Threshold for stamp finding on readout 1 float thresh2, // Threshold for stamp finding on readout 2 float stampSpacing, // Spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window float sysError, // Relative systematic error in images float skyError, // Relative systematic error in images @@ -82,4 +82,11 @@ ) { + PS_ASSERT_PTR_NON_NULL(stamps, false); + PM_ASSERT_READOUT_NON_NULL(ro1, false); + PM_ASSERT_READOUT_NON_NULL(ro2, false); + PS_ASSERT_IMAGE_NON_NULL(subMask, false); + PS_ASSERT_IMAGE_NON_NULL(variance, false); + PS_ASSERT_PTR_NON_NULL(region, false); + psTrace("psModules.imcombine", 3, "Finding stamps...\n"); @@ -87,5 +94,5 @@ *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, - size, footprint, stampSpacing, sysError, skyError, mode); + size, footprint, stampSpacing, normFrac, sysError, skyError, mode); if (!*stamps) { psError(psErrorCodeLast(), false, "Unable to find stamps."); @@ -96,5 +103,5 @@ psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) { + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -110,4 +117,5 @@ const pmReadout *ro1, const pmReadout *ro2, // Input images int stride, // Size for convolution patches + float normFrac, // Fraction of window for normalisation window float sysError, // Systematic error in images float skyError, // Systematic error in images @@ -158,4 +166,8 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false); PS_ASSERT_INT_NONNEGATIVE(stride, false); + if (isfinite(normFrac)) { + PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false); + } if (isfinite(sysError)) { PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false); @@ -210,4 +222,26 @@ } +bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) { + + if (!readout) return true; + + psImage *image = readout->image; + psImage *mask = readout->mask; + psImage *variance = readout->variance; + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue; + bool valid = false; + valid = isfinite(image->data.F32[y][x]); + if (variance) { + valid &= isfinite(variance->data.F32[y][x]); + } + if (valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal; + } + } + + return true; +} bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, @@ -282,5 +316,5 @@ } - if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, kernelError, + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) { psFree(kernels); @@ -289,5 +323,10 @@ } - psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0, + pmSubtractionMaskInvalid(ro1, maskVal); + pmSubtractionMaskInvalid(ro2, maskVal); + + psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels + + psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0, badFrac, mode); // Subtraction mask if (!subMask) { @@ -348,9 +387,9 @@ int inner, int ringsOrder, int binning, float penalty, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, - int iter, float rej, float sysError, float skyError, float kernelError, psImageMaskType maskVal, - psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac, - float badFrac, pmSubtractionMode subMode) + int iter, float rej, float normFrac, float sysError, float skyError, + float kernelError, psImageMaskType maskVal, psImageMaskType maskBad, + psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode) { - if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, skyError, kernelError, + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) { return false; @@ -411,5 +450,5 @@ // Putting important variable declarations here, since they are freed after a "goto" if there is an error. psImage *subMask = NULL; // Mask for subtraction - psRegion *region = NULL; // Iso-kernel region + psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region psString regionString = NULL; // String for region pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF @@ -424,6 +463,10 @@ memCheck("start"); - subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint, - badFrac, subMode); + pmSubtractionMaskInvalid(ro1, maskVal); + pmSubtractionMaskInvalid(ro2, maskVal); + + psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels + + subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode); if (!subMask) { psError(psErrorCodeLast(), false, "Unable to generate subtraction mask."); @@ -435,17 +478,21 @@ // Get region of interest int xRegions = 1, yRegions = 1; // Number of iso-kernel regions - float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions + float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions if (isfinite(regionSize) && regionSize != 0.0) { - xRegions = numCols / regionSize + 1; - yRegions = numRows / regionSize + 1; - xRegionSize = (float)numCols / (float)xRegions; - yRegionSize = (float)numRows / (float)yRegions; - region = psRegionAlloc(NAN, NAN, NAN, NAN); - } - + xRegions = (bounds.x1 - bounds.x0) / regionSize + 1; + yRegions = (bounds.y1 - bounds.y0) / regionSize + 1; + xRegionSize = (float)(bounds.x1 - bounds.x0) / (float)xRegions; + yRegionSize = (float)(bounds.y1 - bounds.y0) / (float)yRegions; + } else { + xRegionSize = bounds.x1 - bounds.x0; + yRegionSize = bounds.y1 - bounds.y0; + } + + // General background subtraction and measurement of stamp threshold float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images { - psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background if (ro1) { + psStatsInit(bg); if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics."); @@ -453,7 +500,9 @@ goto MATCH_ERROR; } - stampThresh1 = bg->robustMedian + threshold * bg->robustStdev; + stampThresh1 = threshold * bg->robustStdev; + psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); } if (ro2) { + psStatsInit(bg); if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics."); @@ -461,8 +510,50 @@ goto MATCH_ERROR; } - stampThresh2 = bg->robustMedian + threshold * bg->robustStdev; - } + stampThresh2 = threshold * bg->robustStdev; + psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); + } psFree(bg); } + + // Just in case the iso-kernel region doesn't cover the entire image... + if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE || + subMode == PM_SUBTRACTION_MODE_DUAL) { + if (!conv1->image) { + conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv1->image, NAN); + if (ro1->variance) { + if (!conv1->variance) { + conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv1->variance, NAN); + } + if (subMask) { + if (!conv1->mask) { + conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } + psImageInit(conv1->mask, maskBad); + } + } + if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_UNSURE || + subMode == PM_SUBTRACTION_MODE_DUAL) { + if (!conv2->image) { + conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv2->image, NAN); + if (ro2->variance) { + if (!conv2->variance) { + conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv2->variance, NAN); + } + if (subMask) { + if (!conv2->mask) { + conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } + psImageInit(conv2->mask, maskBad); + } + } + // Iterate over iso-kernel regions @@ -471,25 +562,27 @@ psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n", j * xRegions + i + 1, xRegions * yRegions); - if (region) { - *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize), - (int)(j * yRegionSize), (int)((j + 1) * yRegionSize)); - psFree(regionString); - regionString = psRegionToString(*region); - psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n", - regionString, numCols, numRows); - } + *region = psRegionSet(bounds.x0 + (int)(i * xRegionSize), + bounds.x0 + (int)((i + 1) * xRegionSize), + bounds.y0 + (int)(j * yRegionSize), + bounds.y0 + (int)((j + 1) * yRegionSize)); + psFree(regionString); + regionString = psRegionToString(*region); + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n", + regionString, numCols, numRows); if (stampsName && strlen(stampsName) > 0) { stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size, - footprint, stampSpacing, sysError, skyError, subMode); + footprint, stampSpacing, normFrac, + sysError, skyError, subMode); } else if (sources) { stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size, - footprint, stampSpacing, sysError, skyError, subMode); + footprint, stampSpacing, normFrac, + sysError, skyError, subMode); } // We get the stamps here; we will also attempt to get stamps at the first iteration, but it // doesn't matter. - if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2, - stampSpacing, sysError, skyError, size, footprint, subMode)) { + if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2, + stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) { goto MATCH_ERROR; } @@ -498,4 +591,5 @@ // generate the window function from the set of stamps if (!pmSubtractionStampsGetWindow(stamps, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get stamp window."); goto MATCH_ERROR; } @@ -503,6 +597,7 @@ // Define kernel basis functions if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold, penalty, subMode); + kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, + optFWHMs, optOrder, stamps, footprint, + optThreshold, penalty, bounds, subMode); if (!kernels) { psErrorClear(); @@ -513,5 +608,5 @@ // Not an ISIS/GUNK kernel, or the optimum kernel search failed kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder, penalty, subMode); + inner, binning, ringsOrder, penalty, bounds, subMode); // pmSubtractionVisualShowKernels(kernels); } @@ -563,6 +658,6 @@ if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, - stampThresh1, stampThresh2, stampSpacing, sysError, skyError, - size, footprint, subMode)) { + stampThresh1, stampThresh2, stampSpacing, normFrac, + sysError, skyError, size, footprint, subMode)) { goto MATCH_ERROR; } @@ -570,4 +665,5 @@ // generate the window function from the set of stamps if (!pmSubtractionStampsGetWindow(stamps, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get stamps window."); goto MATCH_ERROR; } Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.h (revision 26747) @@ -39,4 +39,5 @@ int iter, ///< Rejection iterations float rej, ///< Rejection threshold + float normFrac, ///< Fraction of flux in window for normalisation window float sysError, ///< Relative systematic error in images float skyError, ///< Relative systematic error in images Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionParams.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionParams.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionParams.c (revision 26747) @@ -204,5 +204,6 @@ int spatialOrder, const psVector *fwhms, int maxOrder, const pmSubtractionStampList *stamps, int footprint, - float tolerance, float penalty, pmSubtractionMode mode) + float tolerance, float penalty, psRegion bounds, + pmSubtractionMode mode) { if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { @@ -232,5 +233,5 @@ psVectorInit(orders, maxOrder); pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, - penalty, mode); // Kernels + penalty, bounds, mode); // Kernels psFree(orders); psFree(kernels->description); @@ -483,10 +484,10 @@ if (type == PM_SUBTRACTION_KERNEL_ISIS) { - // XXX in r26035, this code was just wrong. we had: - - // psKernel *subtract = kernels->preCalc->data[0] - - // but, kernels->preCalc was an array of psArray, not an array of kernels. It is now - // an array of pmSubtractionKernelPreCalc. + // XXX in r26035, this code was just wrong. we had: + + // psKernel *subtract = kernels->preCalc->data[0] + + // but, kernels->preCalc was an array of psArray, not an array of kernels. It is now + // an array of pmSubtractionKernelPreCalc. pmSubtractionKernelPreCalc *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionParams.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionParams.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionParams.h (revision 26747) @@ -17,5 +17,6 @@ float tolerance, ///< Maximum difference in chi^2 float penalty, ///< Penalty for wideness - pmSubtractionMode mode // Mode for subtraction + psRegion bounds, ///< Bounds of validity + pmSubtractionMode mode ///< Mode for subtraction ); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.c (revision 26747) @@ -176,5 +176,6 @@ pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, - int footprint, float spacing, float sysErr, float skyErr) + int footprint, float spacing, float normFrac, + float sysErr, float skyErr) { pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return @@ -217,4 +218,6 @@ list->flux = NULL; list->window = NULL; + list->normFrac = normFrac; + list->normWindow = 0; list->footprint = footprint; list->sysErr = sysErr; @@ -239,4 +242,5 @@ out->window = psMemIncrRefCounter(in->window); out->footprint = in->footprint; + out->normWindow = in->normWindow; for (int i = 0; i < num; i++) { @@ -308,4 +312,5 @@ stamp->matrix = NULL; stamp->vector = NULL; + stamp->norm = NAN; return stamp; @@ -316,6 +321,6 @@ const psImage *image2, const psImage *subMask, const psRegion *region, float thresh1, float thresh2, - int size, int footprint, float spacing, float sysErr, float skyErr, - pmSubtractionMode mode) + int size, int footprint, float spacing, float normFrac, + float sysErr, float skyErr, pmSubtractionMode mode) { if (!image1 && !image2) { @@ -373,5 +378,6 @@ if (!stamps) { - stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr); + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, + normFrac, sysErr, skyErr); } @@ -407,6 +413,4 @@ // Take stamps off the top of the (sorted) list for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) { - int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre - // Chop off the top of the list xList->n = j; @@ -414,9 +418,12 @@ fluxList->n = j; +#if 0 // Fish around a bit to see if we can find a pixel that isn't masked + // This is not a good idea if we're using the window feature psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n", i, xCentre, yCentre); // Search bounds + int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre int search = footprint - size; // Search radius int xMin = PS_MAX(border, xCentre - search); @@ -427,8 +434,10 @@ goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, subMask, xMin, xMax, yMin, yMax, numCols, numRows, border); - - // XXX reset for a test: - xStamp = xList->data.F32[j]; - yStamp = yList->data.F32[j]; +#else + // Only search the exact centre pixel + goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, + subMask, xList->data.F32[j], xList->data.F32[j], + yList->data.F32[j], yList->data.F32[j], numCols, numRows, border); +#endif // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp); } @@ -486,5 +495,6 @@ const psImage *image, const psImage *subMask, const psRegion *region, int size, int footprint, - float spacing, float sysErr, float skyErr, pmSubtractionMode mode) + float spacing, float normFrac, float sysErr, float skyErr, + pmSubtractionMode mode) { @@ -507,5 +517,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, region, footprint, spacing, - sysErr, skyErr); // Stamp list + normFrac, sysErr, skyErr); // Stamp list int numStamps = stamps->num; // Number of stamps @@ -612,5 +622,5 @@ PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - int size = kernelSize + stamps->footprint; // Size of postage stamps + int size = stamps->footprint; // Size of postage stamps psFree (stamps->window); @@ -641,13 +651,15 @@ // storage vector for flux data psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); - - double maxValue = 0.0; + psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); + psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); // generate the window pixels + double sum = 0.0; // Sum inside the window + float maxValue = 0.0; // Maximum value, for normalisation for (int y = -size; y <= size; y++) { for (int x = -size; x <= size; x++) { - flux->n = 0; + flux1->n = 0; + flux2->n = 0; for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -656,17 +668,50 @@ if (!stamp->image2) continue; - psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]); - psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]); + psVectorAppend(flux1, stamp->image1->kernel[y][x] / norm1->data.F32[i]); + psVectorAppend(flux2, stamp->image2->kernel[y][x] / norm2->data.F32[i]); } psStatsInit (stats); - if (!psVectorStats (stats, flux, NULL, NULL, 0)) { + if (!psVectorStats (stats, flux1, NULL, NULL, 0)) { psAbort ("failed to generate stats"); } - stamps->window->kernel[y][x] = stats->robustMedian; - if (maxValue < stats->robustMedian) { - maxValue = stats->robustMedian; - } - } + float f1 = stats->robustMedian; + psStatsInit (stats); + if (!psVectorStats (stats, flux2, NULL, NULL, 0)) { + psAbort ("failed to generate stats"); + } + float f2 = stats->robustMedian; + + stamps->window->kernel[y][x] = f1 + f2; + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) { + sum += stamps->window->kernel[y][x]; + } + maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]); + } + } + + psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n", + sum, (1.0 - stamps->normFrac) * sum); + bool done = false; + for (int radius = 1; radius <= size && !done; radius++) { + double within = 0.0; + for (int y = -radius; y <= radius; y++) { + for (int x = -radius; x <= radius; x++) { + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { + within += stamps->window->kernel[y][x]; + } + } + } + psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within); + if (within > (1.0 - stamps->normFrac) * sum) { + stamps->normWindow = radius; + done = true; + } + } + + psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow); + if (stamps->normWindow == 0 || stamps->normWindow >= size) { + psError(PS_ERR_UNKNOWN, true, "Unable to determine normalisation window size."); + return false; } @@ -678,5 +723,5 @@ } -# ifdef TESTING +#if 0 { psFits *fits = psFitsOpen ("window.fits", "w"); @@ -684,8 +729,9 @@ psFitsClose (fits); } -# endif +#endif psFree (stats); - psFree (flux); + psFree (flux1); + psFree(flux2); psFree (norm1); psFree (norm2); @@ -694,5 +740,5 @@ bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, - psImage *variance, int kernelSize) + psImage *variance, int kernelSize, psRegion bounds) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -711,5 +757,4 @@ } - int numCols = image1->numCols, numRows = image1->numRows; // Size of images int size = kernelSize + stamps->footprint; // Size of postage stamps @@ -720,14 +765,10 @@ } - if (isnan(stamp->xNorm)) { - stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols; - } - if (isnan(stamp->yNorm)) { - stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows; - } + p_pmSubtractionPolynomialNormCoords(&stamp->xNorm, &stamp->yNorm, stamp->x, stamp->y, + bounds.x0, bounds.x1, bounds.y0, bounds.y1); int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates - if (x < size || x > numCols - size || y < size || y > numRows - size) { - psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y); + if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) { + psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y); return false; } @@ -788,5 +829,6 @@ pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, float sysErr, float skyErr, + int size, int footprint, float spacing, + float normFrac, float sysErr, float skyErr, pmSubtractionMode mode) { @@ -819,5 +861,6 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, - footprint, spacing, sysErr, skyErr, mode); // Stamps + footprint, spacing, normFrac, + sysErr, skyErr, mode); // Stamps psFree(x); psFree(y); @@ -833,6 +876,6 @@ pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, float sysErr, float skyErr, - pmSubtractionMode mode) + int size, int footprint, float spacing, float normFrac, + float sysErr, float skyErr, pmSubtractionMode mode) { PS_ASSERT_STRING_NON_EMPTY(filename, NULL); @@ -851,5 +894,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint, - spacing, sysErr, skyErr, mode); + spacing, normFrac, sysErr, skyErr, mode); psFree(data); Index: branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.h (revision 26747) @@ -25,4 +25,6 @@ int footprint; ///< Half-size of stamps psKernel *window; ///< window function generated from ensemble of stamps + float normFrac; ///< Fraction of flux in window for normalisation window + int normWindow; ///< Size of window for measuring normalisation float sysErr; ///< Systematic error float skyErr; ///< increase effective readnoise @@ -30,11 +32,13 @@ /// Allocate a list of stamps -pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image - int numRows, // Number of rows in image - const psRegion *region, // Region for stamps, or NULL - int footprint, // Half-size of stamps - float spacing, // Rough average spacing between stamps - float sysErr, // Relative systematic error or NAN - float skyErr // Relative systematic error or NAN +pmSubtractionStampList *pmSubtractionStampListAlloc( + int numCols, // Number of columns in image + int numRows, // Number of rows in image + const psRegion *region, // Region for stamps, or NULL + int footprint, // Half-size of stamps + float spacing, // Rough average spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysErr, // Relative systematic error or NAN + float skyErr // Relative systematic error or NAN ); @@ -78,4 +82,5 @@ psImage *matrix; ///< Least-squares matrix, or NULL psVector *vector; ///< Least-squares vector, or NULL + double norm; ///< Normalisation difference pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp; @@ -85,31 +90,35 @@ /// Find stamps on an image -pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL - const psImage *image1, ///< Image for which to find stamps - const psImage *image2, ///< Image for which to find stamps - const psImage *mask, ///< Mask, or NULL - const psRegion *region, ///< Region to search, or NULL - float thresh1, ///< Threshold for stamps in image 1 - float thresh2, ///< Threshold for stamps in image 2 - int size, ///< Kernel half-size - int footprint, ///< Half-size for stamps - float spacing, ///< Rough spacing for stamps - float sysErr, ///< Relative systematic error in images - float skyErr, ///< Relative systematic error in images - pmSubtractionMode mode ///< Mode for subtraction +pmSubtractionStampList *pmSubtractionStampsFind( + pmSubtractionStampList *stamps, ///< Output stamps, or NULL + const psImage *image1, ///< Image for which to find stamps + const psImage *image2, ///< Image for which to find stamps + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float thresh1, ///< Threshold for stamps in image 1 + float thresh2, ///< Threshold for stamps in image 2 + int size, ///< Kernel half-size + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysErr, ///< Relative systematic error in images + float skyErr, ///< Relative systematic error in images + pmSubtractionMode mode ///< Mode for subtraction ); /// Set stamps based on a list of x,y -pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp - const psVector *y, ///< y coordinates for each stamp - const psImage *image, ///< Image for flux of stamp - const psImage *mask, ///< Mask, or NULL - const psRegion *region, ///< Region to search, or NULL - int size, ///< Kernel half-size - int footprint, ///< Half-size for stamps - float spacing, ///< Rough spacing for stamps - float sysErr, ///< Systematic error in images - float skyErr, ///< Systematic error in images - pmSubtractionMode mode ///< Mode for subtraction +pmSubtractionStampList *pmSubtractionStampsSet( + const psVector *x, ///< x coordinates for each stamp + const psVector *y, ///< y coordinates for each stamp + const psImage *image, ///< Image for flux of stamp + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + int size, ///< Kernel half-size + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysErr, ///< Systematic error in images + float skyErr, ///< Systematic error in images + pmSubtractionMode mode ///< Mode for subtraction ); @@ -123,4 +132,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window float sysErr, ///< Systematic error in images float skyErr, ///< Systematic error in images @@ -137,4 +147,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window float sysErr, ///< Systematic error in images float skyErr, ///< Systematic error in images @@ -142,5 +153,9 @@ ); -bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize); +/// Calculate the window and normalisation window from the stamps +bool pmSubtractionStampsGetWindow( + pmSubtractionStampList *stamps, ///< List of stamps + int kernelSize ///< Half-size of kernel + ); /// Extract stamps from the images @@ -149,5 +164,6 @@ psImage *image2, ///< Input image (or NULL) psImage *variance, ///< Variance map - int kernelSize ///< Kernel half-size + int kernelSize, ///< Kernel half-size + psRegion bounds ///< Bounds of validity ); Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.c (revision 26747) @@ -27,4 +27,7 @@ psFree (detections->oldPeaks); psFree (detections->oldFootprints); + + psFree (detections->newSources); + psFree (detections->allSources); return; } @@ -40,4 +43,6 @@ detections->oldPeaks = NULL; detections->oldFootprints = NULL; + detections->newSources = NULL; + detections->allSources = NULL; detections->last = 0; Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.h (revision 26747) @@ -21,7 +21,9 @@ typedef struct { psArray *footprints; // collection of footprints in the image + psArray *oldFootprints; // collection of footprints previously found psArray *peaks; // collection of all peaks contained by the footprints psArray *oldPeaks; // collection of all peaks previously found - psArray *oldFootprints; // collection of footprints previously found + psArray *newSources; // collection of sources + psArray *allSources; // collection of sources int last; } pmDetections; Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.c (revision 26747) @@ -275,5 +275,5 @@ // psphot-specific function which applies the recipe values // only apply selection to sources within specified region -pmPSFClump pmSourcePSFClump(psRegion *region, psArray *sources, psMetadata *recipe) +pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX) { psTrace("psModules.objects", 10, "---- begin ----\n"); @@ -285,33 +285,7 @@ PS_ASSERT_PTR_NON_NULL(sources, errorClump); - PS_ASSERT_PTR_NON_NULL(recipe, errorClump); - - bool status = true; // Status of MD lookup - float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); - if (!status) { - PSF_SN_LIM = 0; - } - float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); - if (!status) { - PSF_CLUMP_GRID_SCALE = 0.1; - } // find the sigmaX, sigmaY clump { - psF32 SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX"); - if (!status) { - psWarning("MOMENTS_SX_MAX not set in recipe"); - SX_MAX = 10.0; - } - psF32 SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX"); - if (!status) { - psWarning("MOMENTS_SY_MAX not set in recipe"); - SY_MAX = 10.0; - } - psF32 AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); - if (!status) { - psWarning("MOMENTS_AR_MAX not set in recipe"); - AR_MAX = 3.0; - } psF32 AR_MIN = 1.0 / AR_MAX; @@ -399,9 +373,6 @@ psfClump.nSigma = stats->sampleStdev; - const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP"); - if (keep_psf_clump) - { - psMetadataAdd(recipe, PS_LIST_TAIL, - "PSF_CLUMP", PS_DATA_IMAGE, "Image of PSF coefficients", splane); + if (savedImage) { + *savedImage = psMemIncrRefCounter(splane); } psFree (splane); @@ -530,10 +501,9 @@ *****************************************************************************/ -bool pmSourceRoughClass(psRegion *region, psArray *sources, psMetadata *recipe, pmPSFClump clump, psImageMaskType maskSat) +bool pmSourceRoughClass(psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_NSIGMA, pmPSFClump clump, psImageMaskType maskSat) { psTrace("psModules.objects", 10, "---- begin ----"); PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); int Nsat = 0; @@ -548,13 +518,4 @@ psVector *starsn_peaks = psVectorAllocEmpty (sources->n, PS_TYPE_F32); psVector *starsn_moments = psVectorAllocEmpty (sources->n, PS_TYPE_F32); - - // get basic parameters, or set defaults - bool status; - float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM"); - if (!status) PSF_SN_LIM = 20.0; - float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); - if (!status) PSF_CLUMP_NSIGMA = 1.5; - - // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.h (revision 26747) @@ -176,13 +176,15 @@ * * The return value indicates the success (TRUE) of the operation. - * - * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD). - * XXX: Save the clump parameters on the Metadata (TBD) - * - */ + */ + pmPSFClump pmSourcePSFClump( + psImage **savedImage, psRegion *region, ///< restrict measurement to specified region psArray *source, ///< The input pmSource - psMetadata *metadata ///< Contains classification parameters + float PSF_SN_LIM, + float PSF_CLUMP_GRID_SCALE, + psF32 SX_MAX, + psF32 SY_MAX, + psF32 AR_MAX ); @@ -200,5 +202,6 @@ psRegion *region, ///< restrict measurement to specified region psArray *sources, ///< The input pmSources - psMetadata *metadata, ///< Contains classification parameters + float PSF_SN_LIM, ///< min S/N for source to be used for PSF model + float PSF_CLUMP_NSIGMA, ///< size of region around peak of clump for PSF stars pmPSFClump clump, ///< Statistics about the PSF clump psImageMaskType maskSat ///< Mask value for saturated pixels Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmSourceIO.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmSourceIO.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmSourceIO.c (revision 26747) @@ -40,4 +40,5 @@ #include "pmPSF.h" #include "pmModel.h" +#include "pmDetections.h" #include "pmSource.h" #include "pmModelClass.h" @@ -344,10 +345,17 @@ // if sources is NULL, write out an empty table - // input / output sources are stored on the readout->analysis as "PSPHOT.SOURCES" -- a better name might be something like PM_SOURCE_DATA - psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + // input / output sources are stored on the readout->analysis as "PSPHOT.DETECTIONS" + + psArray *sources = NULL; + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (detections) { + sources = detections->allSources; + } if (!sources) { + detections = pmDetectionsAlloc(); sources = psArrayAlloc(0); - psMetadataAddArray(readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "Blank array of sources", sources); - psFree(sources); // Held onto by the metadata, so we can continue to use + detections->allSources = sources; + psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "Blank array of sources", detections); + psFree(detections); // Held onto by the metadata, so we can continue to use } @@ -1031,6 +1039,9 @@ } readout->data_exists = true; - status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "input sources", sources); - psFree (sources); + + pmDetections *detections = pmDetectionsAlloc(); + detections->allSources = sources; + status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "input sources", detections); + psFree (detections); return true; } @@ -1124,9 +1135,10 @@ bool status; - // select the psf of interest - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); - if (!psf) return false; - return true; -} - - + // select the detections of interest + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (!detections) return false; + if (!detections->allSources) return false; + return true; +} + + Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotApResid.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotApResid.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotApResid.c (revision 26747) @@ -34,4 +34,5 @@ #include "pmPSF.h" #include "pmModel.h" +#include "pmDetections.h" #include "pmSource.h" #include "pmSourcePlots.h" @@ -53,4 +54,5 @@ PS_ASSERT_PTR_NON_NULL(layout, false); + bool status; Graphdata graphdata; KapaSection section; @@ -61,7 +63,9 @@ pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); - if (sources == NULL) - return false; + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (detections == NULL) return false; + + psArray *sources = detections->allSources; + if (sources == NULL) return false; int kapa = pmKapaOpen (false); Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotMoments.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotMoments.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotMoments.c (revision 26747) @@ -37,4 +37,5 @@ #include "pmPSF.h" #include "pmModel.h" +#include "pmDetections.h" #include "pmSource.h" #include "pmSourcePlots.h" @@ -54,4 +55,5 @@ PS_ASSERT_PTR_NON_NULL(layout, false); + bool status; Graphdata graphdata; KapaSection section; @@ -62,7 +64,9 @@ pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); - if (sources == NULL) - return false; + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (detections == NULL) return false; + + psArray *sources = detections->allSources; + if (sources == NULL) return false; int kapa = pmKapaOpen (false); Index: branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotPSFModel.c =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotPSFModel.c (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotPSFModel.c (revision 26747) @@ -37,4 +37,5 @@ #include "pmPSF.h" #include "pmModel.h" +#include "pmDetections.h" #include "pmSource.h" #include "pmSourcePlots.h" @@ -56,4 +57,5 @@ PS_ASSERT_PTR_NON_NULL(layout, false); + bool status; Graphdata graphdata; KapaSection section; @@ -64,7 +66,9 @@ pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); - if (sources == NULL) - return false; + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (detections == NULL) return false; + + psArray *sources = detections->allSources; + if (sources == NULL) return false; int kapa = pmKapaOpen (false); Index: branches/eam_branches/psModules.stack.20100120/src/psmodules.h =================================================================== --- branches/eam_branches/psModules.stack.20100120/src/psmodules.h (revision 26685) +++ branches/eam_branches/psModules.stack.20100120/src/psmodules.h (revision 26747) @@ -29,4 +29,5 @@ #include #include +#include #include