Index: /trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 16478) +++ /trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 16479) @@ -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,5 +94,9 @@ 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 16478) +++ /trunk/psModules/src/imcombine/pmStack.c (revision 16479) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-13 02:55:33 $ + * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-14 23:33:09 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -486,8 +486,4 @@ 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 16478) +++ /trunk/psModules/src/imcombine/pmStackReject.c (revision 16479) @@ -11,54 +11,41 @@ #define PIXEL_LIST_BUFFER 100 // Number of pixels to add to list at a time -psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold, - const psArray *subRegions, const psArray *kernels) +psPixels *pmStackReject(const psPixels *in, float threshold, const psArray *regions, + const psArray *solutions, const pmSubtractionKernels *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(subRegions, NULL); - PS_ASSERT_ARRAY_NON_NULL(kernels, NULL); - PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, 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); // Get the original image size - int numRegions = subRegions->n; // Number of regions + int numRegions = regions->n; // Number of regions int numCols = 0, numRows = 0; // Size of original image int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image --- should be 0,0 - int size = 0; // Size of kernel for (int i = 0; i < numRegions; i++) { - psRegion *subRegion = subRegions->data[i]; // Region of interest - if (subRegion->x0 < minCols) { - minCols = subRegion->x0; + psRegion *region = regions->data[i]; // Region of interest + if (region->x0 < minCols) { + minCols = region->x0; } - if (subRegion->y0 < minRows) { - minRows = subRegion->y0; + if (region->y0 < minRows) { + minRows = region->y0; } - if (subRegion->x1 > numCols) { - numCols = subRegion->x1; + if (region->x1 > numCols) { + numCols = region->x1; } - 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 (region->y1 > numRows) { + numRows = region->y1; } } - - // 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); + 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; } - psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1), - 0x01); // Mask + psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve psFree(mask); @@ -68,15 +55,7 @@ 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 - if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 || - region->y0 > valid->y1 || region->y1 < valid->y0)) { - // Region is outside of our sub-image - continue; - } - pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest - if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, true)) { + psRegion *region = regions->data[i]; // Region of interest + if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernels, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); psFree(convRO); @@ -89,10 +68,8 @@ // Image of the kernel at the centre of the region - 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) { + 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) { psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); psFree(convRO); @@ -101,23 +78,13 @@ } float sum = 0.0; - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - sum += image->data.F32[y][x]; + for (int y = 0; y < kernel->numRows; y++) { + for (int x = 0; x < kernel->numCols; x++) { + sum += kernel->data.F32[y][x]; } } - psFree(image); + psFree(kernel); - // 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; - } - } + psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image + psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32)); } psFree(inRO); @@ -127,6 +94,6 @@ // Threshold the convolved image psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked - for (int y = 0; y < convolved->numRows; y++) { - for (int x = 0; x < convolved->numCols; x++) { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { if (convolved->data.F32[y][x] > threshold) { bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x, y); @@ -136,7 +103,7 @@ psFree(convolved); - // Now, grow the mask to include everything that touches a bad pixel in the convolution - mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0xff); - assert(mask->numCols == numCols && mask->numRows == numRows); + // Now, we want to convolve the original pixels properly + mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols, 0, numRows), 0xff); + int size = kernels->size; // Size of kernels for (int i = 0; i < bad->n; i++) { int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest @@ -148,5 +115,4 @@ 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 16478) +++ /trunk/psModules/src/imcombine/pmStackReject.h (revision 16479) @@ -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 *kernels ///< Array of kernel parameters for each region + const psArray *solutions, ///< Array of solution vectors for image + const pmSubtractionKernels *kernels ///< Kernel parameters ); Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 16478) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 16479) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.79 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-14 23:18:34 $ + * @version $Revision: 1.80 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-14 23:33:09 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -50,9 +50,18 @@ // 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++) { - out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][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)); return out; @@ -681,9 +690,19 @@ PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, 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; + 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; + } } @@ -705,5 +724,4 @@ int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions - int x0 = ro1->col0, y0 = ro1->row0; // Image offset psImage *convImage1 = out1->image; // Convolved image @@ -762,9 +780,4 @@ 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 @@ -794,10 +807,8 @@ 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 + y0 + size + 1 - yNormSize/2.0) / - (float)yNormSize; // Normalised coordinate + float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // 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 + x0 + size + 1 - xNormSize/2.0) / - (float)xNormSize; // Normalised coordinate + float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate // Only generate polynomial values every kernel footprint, since we have already assumed Index: /trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 16478) +++ /trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 16479) @@ -14,5 +14,5 @@ -//#define TESTING +#define TESTING ////////////////////////////////////////////////////////////////////////////////////////////////////////////// Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 16478) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 16479) @@ -180,6 +180,4 @@ 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 16478) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 16479) @@ -39,5 +39,4 @@ 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 16478) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 16479) @@ -306,35 +306,6 @@ inner, binning, ringsOrder, mode); } - - // 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); - } + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", + PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); memCheck("kernels"); @@ -470,4 +441,25 @@ 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 16478) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 16479) @@ -8,9 +8,4 @@ #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