Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 17296) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 17297) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.82 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-22 19:50:56 $ + * @version $Revision: 1.83 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-03 01:57:59 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -472,5 +472,4 @@ switch (kernels->mode) { - case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); @@ -638,46 +637,17 @@ const pmSubtractionKernels *kernels, bool doBG, bool useFFT) { - PS_ASSERT_PTR_NON_NULL(out1, false); - PS_ASSERT_PTR_NON_NULL(ro1, false); - PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); - PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false); - if (ro1->weight) { - PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false); - PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); - } + PM_ASSERT_READOUT_NON_NULL(out1, false); + PM_ASSERT_READOUT_NON_NULL(out2, false); + PM_ASSERT_READOUT_NON_NULL(ro1, false); + PM_ASSERT_READOUT_NON_NULL(ro2, false); + PM_ASSERT_READOUT_IMAGE(ro1, false); + PM_ASSERT_READOUT_IMAGE(ro2, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false); PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false); - if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - PS_ASSERT_PTR_NON_NULL(out2, false); - PS_ASSERT_PTR_NON_NULL(ro2, false); - PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); - PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); - if (ro2->weight) { - PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false); - PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false); - } - } if (subMask) { PS_ASSERT_IMAGE_NON_NULL(subMask, false); PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, ro1->image, false); - } - if (out1->image) { - PS_ASSERT_IMAGE_NON_NULL(out1->image, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(out1->image, ro1->image, false); - PS_ASSERT_IMAGE_TYPE(out1->image, PS_TYPE_F32, false); - } - if (out1->mask) { - PS_ASSERT_IMAGE_NON_NULL(out1->mask, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(out1->mask, ro1->image, false); - PS_ASSERT_IMAGE_TYPE(out1->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGE_NON_NULL(subMask, false); - } - if (out1->weight) { - PS_ASSERT_IMAGE_NON_NULL(out1->weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(out1->weight, ro1->image, false); - PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false); } if (region && psRegionIsNaN(*region)) { @@ -688,51 +658,32 @@ } - const pmReadout *source; // Source for image parameters - switch (kernels->mode) { - case PM_SUBTRACTION_MODE_TARGET: - case PM_SUBTRACTION_MODE_1: - case PM_SUBTRACTION_MODE_DUAL: - source = ro1; - break; - case PM_SUBTRACTION_MODE_2: - source = ro2; - break; - default: - psAbort("Unsupported subtraction mode: %x", kernels->mode); - } - psImage *image = source->image, *weight = source->weight; // Image and weight map to convolve - int numCols = image->numCols, numRows = image->numRows; // Image dimensions - int x0 = source->col0, y0 = source->row0; // Image offset - - psImage *convImage1 = out1->image; // Convolved image - if (!convImage1) { - convImage1 = out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - psImage *convMask = NULL; // Convolved mask image - if (subMask) { - if (!out1->mask) { - out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - } - convMask = out1->mask; - psImageInit(convMask, 0); - } - psImage *convWeight1 = NULL; // Convolved weight (variance) image - if (weight) { - if (!out1->weight) { - out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - convWeight1 = out1->weight; - psImageInit(convWeight1, 0.0); - } - - psImage *convImage2 = NULL; - psImage *convWeight2 = NULL; // Convolved products for dual mode - if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + // Inputs + psImage *image1 = ro1->image, *weight1 = ro1->weight; // Image and weight map from input 1 + psImage *image2 = ro2->image, *weight2 = ro2->weight; // Image and weight map from input 2 + int numCols = image1->numCols, numRows = image1->numRows; // Image dimensions + int x0 = image1->col0, y0 = image1->row0; // Image offset + + // Outputs + psImage *convImage1 = NULL, *convWeight1 = NULL; // Convolved image and weight from input 1 + if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + convImage1 = out1->image; + if (!convImage1) { + convImage1 = out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + if (weight1) { + if (!out1->weight) { + out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight1 = out1->weight; + psImageInit(convWeight1, 0.0); + } + } + psImage *convImage2 = NULL, *convWeight2 = NULL; // Convolved image and weight from input 2 + if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { convImage2 = out2->image; if (!convImage2) { convImage2 = out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); } - convWeight2 = NULL; // Convolved weight (variance) image - if (ro2->weight) { + if (weight2) { if (!out2->weight) { out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); @@ -741,10 +692,18 @@ psImageInit(convWeight2, 0.0); } - if (subMask) { - // Copying mask --- they're the same! + } + psImage *convMask = NULL; // Convolved mask image (common to inputs 1 and 2) + if (subMask) { + if (!out1->mask) { + out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + } + convMask = out1->mask; + psImageInit(convMask, 0); + if (out2->mask) { psFree(out2->mask); - out2->mask = psMemIncrRefCounter(convMask); - } - } + } + out2->mask = psMemIncrRefCounter(convMask); + } + int size = kernels->size; // Half-size of kernel @@ -769,5 +728,4 @@ psMaskType maskTarget; // Mark these pixels as bad when propagating the subtractionMask switch (kernels->mode) { - case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: maskSource = PM_SUBTRACTION_MASK_BAD_1; @@ -806,10 +764,13 @@ psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve - convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image, weight, - subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, false); - - if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, ro2->image, ro2->weight, - subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, true); + if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image1, weight1, + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, + false); + } + if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, image2, weight2, + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, + kernels->mode == PM_SUBTRACTION_MODE_DUAL); } @@ -820,13 +781,8 @@ if (subMask->data.PS_TYPE_MASK_DATA[y][x] & maskTarget) { convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; - convImage1->data.F32[y][x] = NAN; - if (weight) { - convWeight1->data.F32[y][x] = NAN; - } } } } } - } } @@ -835,4 +791,22 @@ psFree(polyValues); + // Copy anything that wasn't convolved + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + out2->image = psMemIncrRefCounter(ro2->image); + out2->weight = psMemIncrRefCounter(ro2->weight); + out2->mask = psMemIncrRefCounter(ro2->mask); + break; + case PM_SUBTRACTION_MODE_2: + out1->image = psMemIncrRefCounter(ro1->image); + out1->weight = psMemIncrRefCounter(ro1->weight); + out1->mask = psMemIncrRefCounter(ro1->mask); + break; + case PM_SUBTRACTION_MODE_DUAL: + break; + default: + psAbort("Should never get here."); + } + return true; } Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 17296) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 17297) @@ -521,5 +521,4 @@ bool status; // Status of least-squares matrix/vector calculation switch (kernels->mode) { - case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, @@ -941,5 +940,4 @@ psArray *convolutions; // Convolution postage stamps for each kernel basis function switch (kernels->mode) { - case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: target = stamp->image2; Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 17296) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 17297) @@ -18,5 +18,4 @@ typedef enum { PM_SUBTRACTION_MODE_ERR, // Error in the mode - PM_SUBTRACTION_MODE_TARGET, // Convolve image 1 to match target PSF PM_SUBTRACTION_MODE_1, // Convolve image 1 PM_SUBTRACTION_MODE_2, // Convolve image 2 Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 17296) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 17297) @@ -101,40 +101,13 @@ psMaskType maskBlank, float badFrac, pmSubtractionMode mode) { - PS_ASSERT_PTR_NON_NULL(conv1, false); - if (mode == PM_SUBTRACTION_MODE_DUAL) { - PS_ASSERT_PTR_NON_NULL(conv2, false); - } - PS_ASSERT_PTR_NON_NULL(ro1, false); - PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); - PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false); - if (ro1->mask) { - PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false); - PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false); - } - if (ro1->weight) { - PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false); - PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); - } - if (ro2) { - PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); - PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false); - if (ro2->mask) { - PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false); - PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false); - } - if (ro2->weight) { - PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false); - PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false); - } - } else if (!stampsName && !sources) { - psError(PS_ERR_UNEXPECTED_NULL, true, - "A list of sources is required for convolving to a target PSF."); - return false; - } + PM_ASSERT_READOUT_NON_NULL(conv1, false); + PM_ASSERT_READOUT_NON_NULL(conv2, false); + + PM_ASSERT_READOUT_NON_NULL(ro1, false); + PM_ASSERT_READOUT_NON_NULL(ro2, false); + PM_ASSERT_READOUT_IMAGE(ro1, false); + PM_ASSERT_READOUT_IMAGE(ro2, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false); + PS_ASSERT_INT_NONNEGATIVE(footprint, false); // regionSize can be just about anything (except maybe negative, but it can be NAN) @@ -182,5 +155,5 @@ } - // Reset the output readout, just in case + // Reset the output readouts, just in case if (conv1->image) { psFree(conv1->image); @@ -195,19 +168,28 @@ conv1->weight = NULL; } + if (conv2->image) { + psFree(conv2->image); + conv2->image = NULL; + } + if (conv2->mask) { + psFree(conv2->mask); + conv2->mask = NULL; + } + if (conv2->weight) { + psFree(conv2->weight); + conv2->weight = NULL; + } // Where does our weight map come from? + // Getting the weight exactly right is not necessary --- it's just used for weighting. psImage *weight = NULL; // Weight image to use - if (ro1->weight && ro2 && ro2->weight) { + if (ro1->weight && ro2->weight) { weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight); } else if (ro1->weight) { weight = psMemIncrRefCounter(ro1->weight); - } else if (ro2) { - if (ro2->weight) { - weight = psMemIncrRefCounter(ro2->weight); - } else { - weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); - } + } else if (ro2->weight) { + weight = psMemIncrRefCounter(ro2->weight); } else { - weight = psMemIncrRefCounter(ro1->image); + weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); } @@ -218,4 +200,5 @@ pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF pmSubtractionKernels *kernels = NULL; // Kernel basis functions + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions @@ -272,5 +255,5 @@ } - if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) { + if (mode == PM_SUBTRACTION_MODE_UNSURE) { // Get backgrounds psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background @@ -303,9 +286,4 @@ break; case PM_SUBTRACTION_MODE_2: - if (mode == PM_SUBTRACTION_MODE_TARGET) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Input PSF is larger than target PSF --- can't match image."); - goto MATCH_ERROR; - } psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); break; @@ -334,17 +312,5 @@ // 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; + psRegion *subRegion; // Region over which subtraction was performed if (region) { subRegion = psMemIncrRefCounter(region); @@ -352,8 +318,19 @@ 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) { + + if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) { + psMetadataAddPtr(conv1->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); + 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 (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) { + psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, + PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); + psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction kernels", mode); psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_DATA_REGION | PS_META_DUPLICATE_OK, Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 17296) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 17297) @@ -50,5 +50,5 @@ /// Determine which image to convolve pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, ///< Stamps that have been extracted - int footprint ///< Stamp half-size + float bg1, float bg2 // Background for each image ); Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 17296) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 17297) @@ -86,5 +86,4 @@ switch (mode) { case PM_SUBTRACTION_MODE_1: - case PM_SUBTRACTION_MODE_TARGET: maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1; break; @@ -585,6 +584,5 @@ // Let pmSubtractionStampsSet take care of the rest of the assertions - const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file - psArray *data = psVectorsReadFromFile(filename, format); + psArray *data = psVectorsReadFromFile(filename, "%f %f"); if (!data) { psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename); @@ -592,5 +590,4 @@ } psVector *x = data->data[0], *y = data->data[1]; // Stamp positions - psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions @@ -598,5 +595,5 @@ psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, image, subMask, region, footprint, + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, NULL, image, subMask, region, footprint, spacing, mode); psFree(data);