Index: trunk/ppSub/src/ppSubReadout.c =================================================================== --- trunk/ppSub/src/ppSubReadout.c (revision 14572) +++ trunk/ppSub/src/ppSubReadout.c (revision 14586) @@ -10,31 +10,5 @@ #include "ppSub.h" -//#define TESTING #define WCS_TOLERANCE 0.001 // Tolerance for WCS - -#include -static void memCheck(const char *where) -{ - return; - psMemBlock **leaks = NULL; - int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true); - size_t largestSize = 0; - psMemId largest = 0; - size_t totalSize = 0; - for (int i = 0; i < numLeaks; i++) { - psMemBlock *mb = leaks[i]; - totalSize += mb->userMemorySize; - if (mb->userMemorySize > largestSize) { - largestSize = mb->userMemorySize; - largest = mb->id; - } - } - psFree(leaks); - fprintf(stderr, "%s:\n", where); - fprintf(stderr, " Memory in use: %zd\n", totalSize); - fprintf(stderr, " Largest block: %ld\n", largest); - fprintf(stderr, " sbrk(): %zd\n", (size_t)sbrk(0)); - return; -} bool ppSubReadout(pmConfig *config, const pmFPAview *view) @@ -75,6 +49,7 @@ psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, config->arguments, "MASK.BLANK"), config); // Mask for blank reg. - const char *stampsFile = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps + const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps + // Generate masks if they don't exist int numCols = input->numCols, numRows = input->numRows; // Image dimensions if (!inRO->mask) { @@ -95,245 +70,27 @@ } - memCheck("start"); - - // Mask for subtraction - psImage *subMask = pmSubtractionMask(refRO->mask, inRO->mask, maskBad, size, footprint); - - memCheck("mask"); - -#if 0 - if (!inRO->weight) { - // Need to ensure input image is completely above zero, so we're not weighting by negative numbers - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - if (!psImageBackground(stats, inRO->image, inRO->mask, maskBad, NULL)) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine statistics of input image."); - return false; - } - float diff = stats->robustMedian - 5.0 * stats->robustStdev; - if (diff < 0.0) { - psLogMsg("ppSub", PS_LOG_INFO, "Adding %f to image background to get positive pixels."); - (void)psBinaryOp(in->image, inRO->image, "+", psScalarAlloc(diff, PS_TYPE_F32)); - } - } -#endif - - // Kernel basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, widths, orders, - inner, binning, ringsOrder); - psArray *stamps = NULL; // Stamps for matching PSF - psVector *solution = NULL; // Solution to match PSF - - memCheck("kernels"); - - int xRegions = 1, yRegions = 1; // Number of iso-kernel regions - float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions - psRegion *region = NULL; // Iso-kernel region - 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); + if (!pmSubtractionMatch(outRO, refRO, inRO, footprint, regionSize, spacing, threshold, stampsName, + NAN, type, size, order, widths, orders, inner, ringsOrder, binning, iter, + rej, maskBad, maskBlank)) { + psError(PS_ERR_UNKNOWN, false, "Unable to match images."); + psFree(outRO); + return false; } - psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images - - // Iterate over iso-kernel regions - for (int j = 0; j < yRegions; j++) { - for (int i = 0; i < xRegions; i++) { - psTrace("ppSub", 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)); - psString string = psRegionToString(*region); - psTrace("ppSub", 3, "Iso-kernel region: %s out of %d,%d\n", string, numCols, numRows); - psFree(string); - } - - int numRejected = -1; // Number of rejected stamps in each iteration - for (int k = 0; k < iter && numRejected != 0; k++) { - psTrace("ppSub", 2, "Iteration %d...\n", k); - - psTrace("ppSub", 3, "Finding stamps...\n"); - if (stampsFile) { - iter = 1; // There is no iterating because we use all the stamps we have - psArray *stampData = psVectorsReadFromFile(stampsFile, "%f %f"); // Stamp positions - psVector *xStamp = stampData->data[0]; // x coordinates - psVector *yStamp = stampData->data[1]; // y coordinates - // Correct for IRAF (unit-offset) positions to C (zero-offset) positions - psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - - stamps = pmSubtractionSetStamps(xStamp, yStamp, NULL, subMask, region); - } else { - stamps = pmSubtractionFindStamps(stamps, refRO->image, subMask, region, - threshold, spacing); - } - if (!stamps) { - psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image."); - goto ERROR; - } - - memCheck(" find stamps"); - - psTrace("ppSub", 3, "Extracting stamps...\n"); - if (!pmSubtractionExtractStamps(stamps, refRO->image, inRO->image, inRO->weight, - footprint, kernels)) { - psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); - goto ERROR; - } - - memCheck(" extract stamps"); - - psTrace("ppSub", 3, "Calculating equation...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - - memCheck(" calculate equation"); - - psTrace("ppSub", 3, "Solving equation...\n"); - solution = pmSubtractionSolveEquation(solution, stamps); - if (!solution) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - - memCheck(" solve equation"); - - psTrace("ppSub", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels); - if (numRejected < 0) { - psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); - goto ERROR; - } - psLogMsg("ppSub", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, k); - - memCheck(" reject stamps"); - } - - if (numRejected > 0) { - psTrace("ppSub", 3, "Solving equation...\n"); - solution = pmSubtractionSolveEquation(solution, stamps); - if (!solution) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - } - psFree(stamps); - stamps = NULL; - - memCheck("solution"); - -#ifdef TESTING - psTrace("ppSub", 2, "Generating diagnostics...\n"); - { - // Generate image with convolution kernels - int fullSize = 2 * size + 1 + 1; // Full size of kernel - psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32); - psImageInit(convKernels, NAN); - for (int j = -2; j <= 2; j++) { - for (int i = -2; i <= 2; i++) { - psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0, - (float)j / 2.0); // Image of the kernel - if (!kernel) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); - psFree(convKernels); - goto ERROR; - } - - if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize, - (j + 2) * fullSize, "=") == 0) { - psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image."); - psFree(kernel); - psFree(convKernels); - goto ERROR; - } - psFree(kernel); - } - } - - // XXX What do we do with this image? - - psFits *kernelFile = psFitsOpen("kernel.fits", "w"); - (void)psFitsWriteImage(kernelFile, NULL, convKernels, 0, NULL); - psFitsClose(kernelFile); - - psFree(convKernels); - } - - { - // Generate images of the kernel components - psMetadata *header = psMetadataAlloc(); // Header - for (int i = 0; i < solution->n; i++) { - psString name = NULL; // Header keyword - psStringAppend(&name, "SOLN%04d", i); - psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]); - psFree(name); - } - psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0); - psFits *kernelFile = psFitsOpen("kernels.fits", "w"); - (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL); - psFitsClose(kernelFile); - psFree(kernelImages); - psFree(header); - } -#endif - - memCheck("diag outputs"); - - psTrace("ppSub", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, refRO->image, refRO->weight, - subMask, maskBlank, region, solution, kernels)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image."); - goto ERROR; - } - - psFree(solution); - solution = NULL; - } - } - - psFree(region); - region = NULL; - psFree(subMask); - subMask = NULL; - - if (!pmSubtractionBorder(convImage, convWeight, convMask, kernels, maskBlank)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set border region of convolved image."); - goto ERROR; - } // Add kernel descrption to header - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, - "Subtraction kernel", kernels->description); - psFree(kernels); - - memCheck("convolution"); + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", + psMetadataLookupStr(NULL, outRO->analysis, "SUBTRACTION.KERNEL.DESCRIPTION")); // Do the subtraction if (reverse) { - outRO->image = (psImage*)psBinaryOp(NULL, convImage, "-", inRO->image); + psBinaryOp(outRO->image, outRO->image, "-", inRO->image); } else { - outRO->image = (psImage*)psBinaryOp(NULL, inRO->image, "-", convImage); + psBinaryOp(outRO->image, inRO->image, "-", outRO->image); } - outRO->mask = (psImage*)psBinaryOp(NULL, convMask, "|", inRO->mask); - if (convWeight) { - outRO->weight = (psImage*)psBinaryOp(NULL, convWeight, "+", inRO->weight); + psBinaryOp(outRO->mask, outRO->mask, "|", inRO->mask); + if (outRO->weight) { + psBinaryOp(outRO->weight, outRO->weight, "+", inRO->weight); } - - memCheck("subtraction"); - - psFree(convImage); - psFree(convMask); - psFree(convWeight); - - memCheck("subtraction freed"); - - outRO->data_exists = true; - outCell->data_exists = true; - outCell->parent->data_exists = true; #ifdef TESTING @@ -349,4 +106,5 @@ if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) { psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output."); + psFree(outRO); return false; } @@ -387,17 +145,6 @@ "Subtraction input", inFile->filename); - memCheck("header"); - psFree(outRO); return true; - -ERROR: - psFree(region); - psFree(subMask); - psFree(kernels); - psFree(stamps); - psFree(solution); - psFree(outRO); - return false; }