Index: trunk/ppSub/src/ppSubArguments.c =================================================================== --- trunk/ppSub/src/ppSubArguments.c (revision 14404) +++ trunk/ppSub/src/ppSubArguments.c (revision 14457) @@ -179,5 +179,5 @@ psMetadataAddStr(arguments, PS_LIST_TAIL, "-refweight", 0, "Referenceweight image", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL); - + psMetadataAddF32(arguments, PS_LIST_TAIL, "-region", 0, "Size of iso-kernel region", NAN); psMetadataAddS32(arguments, PS_LIST_TAIL, "-size", 0, "Kernel half-size (pixels)", 0); psMetadataAddS32(arguments, PS_LIST_TAIL, "-order", 0, "Spatial polynomial order", -1); @@ -232,4 +232,5 @@ } + VALUE_ARG_RECIPE_FLOAT("-region", "REGION.SIZE", F32); VALUE_ARG_RECIPE_INT("-size", "KERNEL.SIZE", S32, 0); VALUE_ARG_RECIPE_INT("-order", "SPATIAL.ORDER", S32, -1); Index: trunk/ppSub/src/ppSubReadout.c =================================================================== --- trunk/ppSub/src/ppSubReadout.c (revision 14404) +++ trunk/ppSub/src/ppSubReadout.c (revision 14457) @@ -10,5 +10,5 @@ #include "ppSub.h" -#define TESTING +//#define TESTING #define WCS_TOLERANCE 0.001 // Tolerance for WCS @@ -52,4 +52,5 @@ int size = psMetadataLookupS32(NULL, config->arguments, "KERNEL.SIZE"); // Kernel half-size int order = psMetadataLookupS32(NULL, config->arguments, "SPATIAL.ORDER"); // Spatial polynomial order + float regionSize = psMetadataLookupF32(NULL, config->arguments, "REGION.SIZE"); // Size of iso-kernel regs float spacing = psMetadataLookupF32(NULL, config->arguments, "STAMP.SPACING"); // Typical stamp spacing int footprint = psMetadataLookupS32(NULL, config->arguments, "STAMP.FOOTPRINT"); // Stamp half-size @@ -119,120 +120,165 @@ memCheck("kernels"); - int numRejected = -1; // Number of rejected stamps in each iteration - for (int i = 0; i < iter && numRejected != 0; i++) { - stamps = pmSubtractionFindStamps(stamps, refRO->image, subMask, threshold, spacing); - if (!stamps) { - psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image."); - goto ERROR; - } - - memCheck(" find stamps"); - - if (!pmSubtractionCalculateEquation(stamps, refRO->image, inRO->image, inRO->weight, - kernels, footprint)) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - - memCheck(" calculate equation"); - - solution = pmSubtractionSolveEquation(solution, stamps); - if (!solution) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - - memCheck(" solve equation"); - - numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, 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, i); - - memCheck(" reject stamps"); - } - - if (numRejected > 0) { - 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"); + 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)) { + xRegions = numCols / regionSize + 1; + yRegions = numRows / regionSize + 1; + xRegionSize = (float)numCols / (float)xRegions; + yRegionSize = (float)numRows / (float)yRegions; + region = psRegionAlloc(NAN, NAN, NAN, NAN); + } + + 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"); + 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, "Calculating equation...\n"); + if (!pmSubtractionCalculateEquation(stamps, refRO->image, inRO->image, inRO->weight, + 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, refRO->image, inRO->image, 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 - { - // 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); - } + 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"); - - psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images - if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, refRO->image, refRO->weight, subMask, - maskBlank, solution, kernels)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image."); - goto ERROR; - } + 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; + } + psFree(kernels); - psFree(solution); memCheck("convolution"); @@ -265,5 +311,5 @@ for (int x = 0; x < outRO->image->numCols; x++) { if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.U8[y][x] & maskBlank)) { - printf("%d %d --> %d\n", x, y, outRO->mask->data.U8[y][x]); + printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.U8[y][x]); } } @@ -315,4 +361,5 @@ ERROR: + psFree(region); psFree(subMask); psFree(kernels);