Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 16604) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 16607) @@ -31,5 +31,5 @@ -#define TESTING // Enable test output +//#define TESTING // Enable test output #define PEAK_FLUX 1.0e3 // Peak flux for each source #define SKY_VALUE 0.0e0 // Sky value for fake image @@ -94,9 +94,5 @@ fake->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE); fake->type = PM_SOURCE_TYPE_STAR; - -#ifdef TESTING - // Required to generate model image fake->psfMag = -2.5 * log10(PEAK_FLUX); -#endif psTrace("psModules.imcombine", 5, "Source %d: %.2f,%.2f\n", Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 16604) +++ trunk/psModules/src/imcombine/pmStack.c (revision 16607) @@ -8,4 +8,6 @@ * @author GLG, MHPCC * + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-22 19:50:56 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -510,4 +512,8 @@ for (int j = 0; j < pixels->n; j++) { int x = pixels->data[j].x, y = pixels->data[j].y; // Coordinates of interest + if (x < 0 || x >= numCols || y < 0 || y >= numRows) { + psWarning("Bad pixel coordinate: %d,%d --- ignored.", x, y); + continue; + } psArray *columns = map->data[y]; // The columns for that row psVector *images = columns->data[x]; // The images for that column Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 16604) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 16607) @@ -11,43 +11,55 @@ #define PIXEL_LIST_BUFFER 100 // Number of pixels to add to list at a time -psPixels *pmStackReject(const psPixels *in, float threshold, const psArray *regions, - const psArray *solutions, const pmSubtractionKernels *kernels) +psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold, + const psArray *subRegions, const psArray *kernels) { PS_ASSERT_PIXELS_NON_NULL(in, NULL); PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(threshold, 0.0, NULL); PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL); - PS_ASSERT_ARRAY_NON_NULL(regions, NULL); - PS_ASSERT_ARRAY_NON_NULL(solutions, NULL); - PS_ASSERT_ARRAYS_SIZE_EQUAL(regions, solutions, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL); + PS_ASSERT_ARRAY_NON_NULL(kernels, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL); // Get the original image size - int numRegions = regions->n; // Number of regions + int numRegions = subRegions->n; // Number of regions int numCols = 0, numRows = 0; // Size of original image int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image int size = 0; // Size of kernel for (int i = 0; i < numRegions; i++) { - psRegion *region = regions->data[i]; // Region of interest - if (region->x0 < minCols) { - minCols = region->x0; + psRegion *subRegion = subRegions->data[i]; // Region of interest + if (subRegion->x0 < minCols) { + minCols = subRegion->x0; } - if (region->y0 < minRows) { - minRows = region->y0; + if (subRegion->y0 < minRows) { + minRows = subRegion->y0; } - if (region->x1 > numCols) { - numCols = region->x1; + if (subRegion->x1 > numCols) { + numCols = subRegion->x1; } - if (region->y1 > numRows) { - numRows = region->y1; + if (subRegion->y1 > numRows) { + numRows = subRegion->y1; + } + + pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest + if (size == 0) { + size = kernel->size; + } else if (kernel->size != size) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d", + size, kernel->size); + return NULL; } } - if (minCols != 0 || minRows != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Some error with image regions --- minimum coordinate is not 0,0"); - return NULL; + + // Adjust the size for the size of the subimage + if (valid) { + minCols = PS_MAX(valid->x0, minCols); + minRows = PS_MAX(valid->y0, minRows); + numCols = PS_MIN(valid->x1, numCols); + numRows = PS_MIN(valid->y1, numRows); } psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows); - psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image + psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1), + 0x01); // Mask psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve psFree(mask); @@ -57,4 +69,6 @@ pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image inRO->image = image; + inRO->col0 = minCols; + inRO->row0 = minRows; for (int i = 0; i < numRegions; i++) { psRegion *region = subRegions->data[i]; // Region of interest @@ -76,8 +90,10 @@ // Image of the kernel at the centre of the region - float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - numCols/2.0) / (float)numCols; - float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - numRows/2.0) / (float)numRows; - psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false); - if (!kernel) { + float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) / + (float)kernel->numCols; + float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) / + (float)kernel->numRows; + psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false); + if (!image) { psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); psFree(convRO); @@ -86,13 +102,23 @@ } float sum = 0.0; - for (int y = 0; y < kernel->numRows; y++) { - for (int x = 0; x < kernel->numCols; x++) { - sum += kernel->data.F32[y][x]; + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + sum += image->data.F32[y][x]; } } - psFree(kernel); + psFree(image); - psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image - psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32)); + // Range for normalisation + int yMin = PS_MAX(minRows, region->y0) - inRO->row0; + int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0; + int xMin = PS_MAX(minCols, region->x0) - inRO->col0; + int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0; + psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n", + sum, xMin, xMax, yMin, yMax); + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + convRO->image->data.F32[y][x] /= sum; + } + } } psFree(inRO); @@ -124,4 +150,5 @@ for (int y = yMin; y <= yMax; y++) { for (int x = xMin; x <= xMax; x++) { + assert(x < mask->numCols && y < mask->numRows); mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff; } Index: trunk/psModules/src/imcombine/pmStackReject.h =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.h (revision 16604) +++ trunk/psModules/src/imcombine/pmStackReject.h (revision 16607) @@ -10,8 +10,8 @@ /// We apply a matched filter to the corresponding mask image, and threshold to find the original pixels psPixels *pmStackReject(const psPixels *in, ///< List of pixels in the convolved image + const psRegion *valid, ///< Valid region to consider float threshold, ///< Threshold on convolved image, 0..1 const psArray *regions, ///< Array of image regions for image - const psArray *solutions, ///< Array of solution vectors for image - const pmSubtractionKernels *kernels ///< Kernel parameters + const psArray *kernels ///< Array of kernel parameters for each region ); Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 16604) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 16607) @@ -4,4 +4,7 @@ * @author GLG, MHPCC * + * @version $Revision: 1.82 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-22 19:50:56 $ + * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -47,18 +50,9 @@ // Take the square of the normal kernel - double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels for (int v = yMin; v <= yMax; v++) { for (int u = xMin; u <= xMax; u++) { - float value = normalKernel->kernel[v][u]; // Value of interest - float value2 = PS_SQR(value); // Value squared - sumNormal += value; - sumVariance += value2; - out->kernel[v][u] = value2; - } - } - - // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel - // This is required to keep the relative scaling between the image and the weight map - psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32)); + out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]); + } + } return out; @@ -687,19 +681,9 @@ PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false); } - if (region) { - if (psRegionIsNaN(*region)) { - psString string = psRegionToString(*region); - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string); - psFree(string); - return false; - } - if (region->x0 < 0 || region->x1 > ro1->image->numCols || - region->y0 < 0 || region->y1 > ro1->image->numRows) { - psString string = psRegionToString(*region); - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)", - string, ro1->image->numCols, ro1->image->numRows); - psFree(string); - return false; - } + if (region && psRegionIsNaN(*region)) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string); + psFree(string); + return false; } @@ -776,4 +760,9 @@ yMax = PS_MIN(region->y1, yMax); } + + // 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); psMaskType maskSource; // Mask these pixels when convolving @@ -803,8 +792,10 @@ for (int j = yMin; j < yMax; j += fullSize) { int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest - float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // Normalised coordinate + float yNorm = 2.0 * (float)(j + y0 + size + 1 - yNormSize/2.0) / + (float)yNormSize; // Normalised coordinate for (int i = xMin; i < xMax; i += fullSize) { int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest - float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate + float xNorm = 2.0 * (float)(i + x0 + size + 1 - xNormSize/2.0) / + (float)xNormSize; // Normalised coordinate // Only generate polynomial values every kernel footprint, since we have already assumed Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 16604) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 16607) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-22 19:24:42 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-22 19:50:56 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 16604) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 16607) @@ -14,5 +14,5 @@ -#define TESTING +//#define TESTING ////////////////////////////////////////////////////////////////////////////////////////////////////////////// Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 16604) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 16607) @@ -180,4 +180,6 @@ kernels->bgOrder = 0; kernels->mode = mode; + kernels->numCols = 0; + kernels->numRows = 0; kernels->solution1 = NULL; kernels->solution2 = NULL; Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 16604) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 16607) @@ -39,4 +39,5 @@ 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 } pmSubtractionKernels; Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 16604) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 16607) @@ -306,6 +306,35 @@ inner, binning, ringsOrder, mode); } - psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", - PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); + + // Add analysis metadata + { + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, + PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); + if (conv2) { + psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, + PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); + } + psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction kernels", mode); + if (conv2) { + psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction kernels", mode); + } + psRegion *subRegion; + if (region) { + subRegion = psMemIncrRefCounter(region); + } else { + subRegion = psRegionAlloc(0, numCols, 0, numRows); + } + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, + PS_DATA_REGION | PS_META_DUPLICATE_OK, + "Region over which subtraction was performed", subRegion); + if (conv2) { + psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, + PS_DATA_REGION | PS_META_DUPLICATE_OK, + "Region over which subtraction was performed", subRegion); + } + psFree(subRegion); + } memCheck("kernels"); @@ -442,25 +471,4 @@ kernels = NULL; -#if 0 - // Put the solution on the metadata - { - psString comment = NULL; // Comment for metadata - psStringAppend(&comment, "Subtraction solution for region %s", regionString); - psMetadataAddVector(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION", - PS_META_DUPLICATE_OK, comment, solution); - psFree(comment); - if (region) { - psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION", - PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region); - } else { - region = psRegionAlloc(0, numCols, 0, numRows); - psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION", - PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region); - psFree(region); - region = NULL; - } - } -#endif - // There is data in the readout now conv1->data_exists = true; Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 16604) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 16607) @@ -8,4 +8,9 @@ #include #include + +#define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Name of kernel in the analysis metadata +#define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Name of subtraction mode in the analysis metadata +#define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Name of subtraction region in the analysis MD + /// Match two images