Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 14019) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 14106) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-06-09 01:06:48 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-07-10 23:54:26 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -24,4 +24,6 @@ #include "pmSubtraction.h" + +//#define CAN_PASS_STATIC_INLINE // Can I pass the address of a static inline function? ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -53,4 +55,5 @@ } +#ifndef CAN_PASS_STATIC_INLINE // Generate an image of the solved kernel // Meant to replace the following function declaration: @@ -61,5 +64,5 @@ // double (*weightFunc)(double value) // Function for weighting // ); -// because even if the weightFunc can't be "static inline", they still appear in the assembly. +// because even though the weightFunc is "static inline", it still appears in the assembly. // TARGET: Kernel, to return (psKernel*) // SOLUTION: Solution to the least-squares problem (psVector*) @@ -75,5 +78,7 @@ \ /* Ensure the subIndex for POIS kernels is what is expected */ \ - assert((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS || \ + assert(((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS && \ + (KERNELS)->type != PM_SUBTRACTION_KERNEL_SPAM && \ + (KERNELS)->type != PM_SUBTRACTION_KERNEL_FRIES) || \ ((KERNELS)->u->data.S32[(KERNELS)->subIndex] == 0 && \ (KERNELS)->v->data.S32[(KERNELS)->subIndex] == 0 && \ @@ -103,4 +108,22 @@ break; \ } \ + /* SPAM and FRIES use the same method */ \ + case PM_SUBTRACTION_KERNEL_SPAM: \ + case PM_SUBTRACTION_KERNEL_FRIES: { \ + int uStart = (KERNELS)->u->data.S32[i]; \ + int uStop = (KERNELS)->uStop->data.S32[i]; \ + int vStart = (KERNELS)->v->data.S32[i]; \ + int vStop = (KERNELS)->vStop->data.S32[i]; \ + for (int v = vStart; v <= vStop; v++) { \ + for (int u = uStart; u <= uStop; u++) { \ + (TARGET)->kernel[v][u] += value; \ + } \ + } \ + if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \ + /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \ + (TARGET)->kernel[0][0] += subValue; \ + } \ + break; \ + } \ case PM_SUBTRACTION_KERNEL_ISIS: { \ psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \ @@ -122,6 +145,5 @@ } \ } - -#if 0 +#else // Generate an image of the solved kernel static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return @@ -139,5 +161,7 @@ // Ensure the subIndex for POIS kernels is what is expected - assert(kernels->type != PM_SUBTRACTION_KERNEL_POIS || + assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS && + kernels->type != PM_SUBTRACTION_KERNEL_SPAM && + kernels->type != PM_SUBTRACTION_KERNEL_FRIES) || (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0 && kernels->xOrder->data.S32[kernels->subIndex] == 0 && @@ -167,4 +191,21 @@ break; } + /* SPAM and FRIES use the same method */ + case PM_SUBTRACTION_KERNEL_SPAM: + case PM_SUBTRACTION_KERNEL_FRIES: + int uStart = kernels->u->data.S32[i]; + int uStop = kernels->uStop->data.S32[i]; + int vStart = kernels->v->data.S32[i]; + int vStop = kernels->vStop->data.S32[i]; + for (int v = vStart; v <= vStop; v++) { + for (int u = uStart; u <= uStop; u++) { + kernel->kernel[v][u] += value; + } + } + if (kernels->spatialOrder > 0 && i != kernels->subIndex) { + /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ + kernel->kernel[0][0] += subValue; + } + break; case PM_SUBTRACTION_KERNEL_ISIS: { psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values @@ -190,4 +231,5 @@ #endif +#ifndef CAN_PASS_STATIC_INLINE // Generate the convolved pixel value // Meant to replace the following function declaration: @@ -199,5 +241,5 @@ // double (*weightFunc)(double value) // Function for weighting // ); -// because even if the weightFunc can't be "static inline", they still appear in the assembly. +// because even though weightFunc is "static inline", it still appears in the assembly. // // TARGET: The value to 'return' (double) @@ -225,5 +267,23 @@ break; \ } \ - case PM_SUBTRACTION_KERNEL_ISIS: { \ + /* Method for SPAM and FRIES is the same */ \ + case PM_SUBTRACTION_KERNEL_SPAM: \ + case PM_SUBTRACTION_KERNEL_FRIES: { \ + int uStart = (KERNELS)->u->data.S32[(INDEX)]; \ + int uStop = (KERNELS)->uStop->data.S32[(INDEX)]; \ + int vStart = (KERNELS)->v->data.S32[(INDEX)]; \ + int vStop = (KERNELS)->vStop->data.S32[(INDEX)]; \ + for (int v = vStart; v <= vStop; v++) { \ + for (int u = uStart; u <= uStop; u++) { \ + (TARGET) += FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \ + } \ + } \ + if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \ + /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \ + (TARGET) += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \ + } \ + break; \ + } \ + case PM_SUBTRACTION_KERNEL_ISIS: { \ psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \ psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \ @@ -231,10 +291,10 @@ double sum = 0.0; /* Accumulated sum from convolution */ \ double sub = 0.0; /* Accumulated sum to subtract */ \ - for (int v = -size; v <= size; v++) { \ - for (int u = -size; u <= size; u++) { \ - sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \ + for (int v = -size, y = (Y) - size; v <= size; v++, y++) { \ + for (int u = -size, x = (X) - size; u <= size; u++, x++) { \ + sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[y][x]; \ /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \ if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \ - sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \ + sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[y][x]; \ } \ } \ @@ -247,6 +307,5 @@ } \ } - -#if 0 +#else // Generate the convolved pixel value static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions @@ -273,4 +332,22 @@ } return value; + } + /* Method for SPAM and FRIES is the same */ + case PM_SUBTRACTION_KERNEL_SPAM: + case PM_SUBTRACTION_KERNEL_FRIES: { + int uStart = kernels->u->data.S32[index]; + int uStop = kernels->uStop->data.S32[index]; + int vStart = kernels->v->data.S32[index]; + int vStop = kernels->vStop->data.S32[index]; + for (int v = vStart; v <= vStop; v++) { + for (int u = uStart; u <= uStop; u++) { + value += weightFunc(polyValue) * image->data.F32[y + v][x + u]; + } + } + if (kernels->spatialOrder > 0 && index != kernels->subIndex) { + /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ + value += weightFunc(-1.0) * image->data.F32[y][x]; + } + break; } case PM_SUBTRACTION_KERNEL_ISIS: { @@ -333,107 +410,4 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -psImage *pmSubtractionMask(const psImage *inMask, const psImage *refMask, psMaskType maskVal, - int size, int footprint) -{ - PS_ASSERT_IMAGE_NON_NULL(inMask, NULL); - PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL); - PS_ASSERT_IMAGE_NON_NULL(refMask, NULL); - PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL); - PS_ASSERT_INT_NONNEGATIVE(size, NULL); - PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); - - // Size of the images - int numCols = inMask->numCols; - int numRows = inMask->numRows; - - // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask - psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask - psImageInit(mask, 0); - - // Dereference for convenience - psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; - psMaskType **inData = inMask->data.PS_TYPE_MASK_DATA; - psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA; - - // Block out a border around the edge of the image - - // Bottom stripe - for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { - for (int x = 0; x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - // Either side - for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { - for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - // Top stripe - for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - - // Mask the bad pixels, and around them - for (int y = footprint + size; y < numRows - footprint - size; y++) { - for (int x = 0; x < numCols; x++) { - if (inData[y][x] & maskVal) { - maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT; - // Block out the entire stamp footprint around this pixel - for (int v = y - footprint; v <= y + footprint; v++) { - for (int u = x - footprint; u <= x + footprint; u++) { - maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; - } - } - } - if (refData[y][x] & maskVal) { - maskData[y][x] |= PM_SUBTRACTION_MASK_REF; - - // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with - // this bad pixel (within 'size'). Then we want to block out with the FOOTPRINT mask - // everything within a footprint's distance of those (within 'footprint'). - - // Bottom stripe - for (int v = PS_MAX(y - footprint - size, 0); v < y - size; v++) { - for (int u = PS_MAX(x - footprint - size, 0); - u <= PS_MIN(x + footprint + size, numCols - 1); u++) { - maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; - } - } - // Middle - for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) { - // Left side - for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) { - maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; - } - // Centre - for (int u = PS_MAX(x - size, 0); u <= PS_MIN(x + size, numCols - 1); u++) { - maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_CONVOLVE; - } - // Right side - for (int u = x + size + 1; u <= PS_MIN(x + footprint + size, numCols - 1); u++) { - maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; - } - } - // Top stripe - for (int v = y + size + 1; v <= PS_MIN(y + footprint + size, numRows - 1); v++) { - for (int u = PS_MAX(x - footprint - size, 0); - u <= PS_MIN(x + footprint + size, numCols - 1); u++) { - maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; - } - } - } - } - } - - return mask; -} - bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input, const psImage *weight, const pmSubtractionKernels *kernels, int footprint) @@ -504,5 +478,5 @@ // Generate the convolutions for (int i = 0; i < numKernels; i++) { -#if 0 +#ifdef CAN_PASS_STATIC_INLINE convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues, imageWeighting); @@ -533,30 +507,13 @@ psFree(polyValues); - // Fill in lower diagonal of symmetric matrix, while checking for bad values - bool bad = false; // Are there bad values? + // Fill in lower diagonal of symmetric matrix for (int i = 0; i < numKernels; i++) { for (int j = 0; j < i; j++) { - if (!isfinite(stampMatrix->data.F64[j][i])) { - bad = true; - } stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i]; } - if (!isfinite(stampMatrix->data.F64[i][bgIndex]) || - !isfinite(stampMatrix->data.F64[i][i]) || - !isfinite(stampVector->data.F64[i])) { - bad = true; - } stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex]; } - if (!isfinite(stampVector->data.F64[bgIndex])) { - bad = true; - } - - if (bad) { - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - psTrace("psModules.imcombine", 3, "Rejecting stamp %d because of bad equation\n", i); - } else { - stamp->status = PM_SUBTRACTION_STAMP_USED; - } + + stamp->status = PM_SUBTRACTION_STAMP_USED; if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { @@ -653,6 +610,6 @@ int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage, - psImage *subMask, const psVector *solution, int footprint, float sigmaRej, - const pmSubtractionKernels *kernels) + psImage *mask, psMaskType badStampMaskVal, const psVector *solution, + int footprint, float sigmaRej, const pmSubtractionKernels *kernels) { PS_ASSERT_ARRAY_NON_NULL(stamps, -1); @@ -662,7 +619,7 @@ PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1); PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -1); - PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); - PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); - PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, subMask, -1); + PS_ASSERT_IMAGE_NON_EMPTY(mask, -1); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, -1); + PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, -1); PS_ASSERT_VECTOR_NON_NULL(solution, -1); PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1); @@ -707,5 +664,5 @@ for (int x = footprint; x <= footprint; x++) { for (int k = 0; k < numKernels; k++) { -#if 0 +#ifdef CAN_PASS_STATIC_INLINE convolvedStamp->data.F32[y + footprint][x + footprint] += convolvePixel(kernels, k, xStamp + x, yStamp + y, refImage, polyValues, @@ -753,9 +710,9 @@ if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) { // Mask out the stamp in the image so you it's not found again - psTrace("psModules.imcombine", 3, "Clipping stamp %d\n", i); + psTrace("psModules.imcombine", 3, "Rejecting stamp %d\n", i); numRejected++; for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { - subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; + mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal; } } @@ -794,5 +751,5 @@ // The appropriate kernel -#if 0 +#ifdef CAN_PASS_STATIC_INLINE psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting); #else @@ -810,13 +767,14 @@ bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, - const psImage *inImage, const psImage *inWeight, const psImage *subMask, - psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels) + const psImage *inImage, const psImage *inWeight, const psImage *inMask, + psMaskType maskVal, psMaskType blank, + const psVector *solution, const pmSubtractionKernels *kernels) { PS_ASSERT_IMAGE_NON_NULL(inImage, false); PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, 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, inImage, false); + if (inMask) { + PS_ASSERT_IMAGE_NON_NULL(inMask, false); + PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, inImage, false); } if (inWeight) { @@ -835,5 +793,5 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false); PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGE_NON_NULL(subMask, false); + PS_ASSERT_IMAGE_NON_NULL(inMask, false); } if (outWeight && *outWeight) { @@ -865,5 +823,5 @@ } psImage *convMask = NULL; // Convolved mask image - if (outMask && subMask) { + if (outMask && inMask) { if (!*outMask) { *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); @@ -896,5 +854,5 @@ 2.0 * (float)(i + size - numCols/2.0) / (float)numCols, 2.0 * (float)(j + size - numRows/2.0) / (float)numRows); -#if 0 +#ifdef CAN_PASS_STATIC_INLINE kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting); if (inWeight) { @@ -914,12 +872,18 @@ // Check and propagate the kernel footprint, if required - if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] & - (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) { - convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; - convImage->data.F32[y][x] = NAN; - if (inWeight) { - convWeight->data.F32[y][x] = NAN; + if (inMask) { + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + convMask->data.PS_TYPE_MASK_DATA[y][x] |= + inMask->data.PS_TYPE_MASK_DATA[y + v][x + u]; + } } - continue; + if (convMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) { + convImage->data.F32[y][x] = NAN; + if (inWeight) { + convWeight->data.F32[y][x] = NAN; + } + continue; + } } Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 14019) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 14106) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-06-09 01:04:02 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-07-10 23:54:26 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -19,22 +19,4 @@ /// @addtogroup imcombine Image Combinations /// @{ - -typedef enum { - PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking - PM_SUBTRACTION_MASK_REF = 0x01, // Reference image is bad - PM_SUBTRACTION_MASK_INPUT = 0x02, // Input image is bad - PM_SUBTRACTION_MASK_CONVOLVE = 0x04, // If convolved, would be bad - PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint - PM_SUBTRACTION_MASK_BORDER = 0x10, // Image border - PM_SUBTRACTION_MASK_REJ = 0x20, // Previously tried as a stamp, and rejected -} pmSubtractionMasks; - -/// Generate a mask for use in the subtraction process -psImage *pmSubtractionMask(const psImage *inMask, ///< Mask for the input image - const psImage *refMask, ///< Mask for the reference image (will be convolved) - psMaskType maskVal, ///< Value to mask out - int size, ///< Half-size of the kernel (pmSubtractionKernels.size) - int footprint ///< Half-size of the kernel footprint - ); /// Calculate the least-squares equation to match the image quality @@ -56,5 +38,6 @@ const psImage *refImage, ///< Reference image psImage *inImage, ///< Input image - psImage *subMask, ///< Subtraction mask + psImage *mask, ///< Mask image + psMaskType badStampMaskVal, ///< Value to use in mask for bad stamp const psVector *solution, ///< Solution vector int footprint, ///< Region to mask if stamp is bad @@ -75,5 +58,6 @@ const psImage *inImage, ///< Input image const psImage *inWeight, ///< Input weight map (or NULL) - const psImage *subMask, ///< Subtraction mask (or NULL) + const psImage *inMask, ///< Input mask (or NULL) + psMaskType maskVal, ///< Value to mask psMaskType blank, ///< Mask value for blank regions const psVector *solution, ///< The solution vector Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14019) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14106) @@ -15,4 +15,6 @@ psFree(kernels->v); psFree(kernels->sigma); + psFree(kernels->uStop); + psFree(kernels->vStop); psFree(kernels->xOrder); psFree(kernels->yOrder); @@ -48,7 +50,9 @@ kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->sigma = NULL; + kernels->uStop = NULL; + kernels->vStop = NULL; kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); - kernels->sigma = NULL; kernels->subIndex = 0; kernels->preCalc = NULL; @@ -187,6 +191,208 @@ } +/// Generate SPAM kernels +pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner, ///< Inner radius to preserve unbinned + int binning ///< Kernel binning factor + ) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + PS_ASSERT_INT_POSITIVE(binning, NULL); + + // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region + int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region + int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region + assert(numOuter * binning + numInner == size); + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1) * + (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions + + psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, + size, spatialOrder); // The kernels + + kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + + psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element + psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element + locations->data.S32[numTotal] = 0; + widths->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + locations->data.S32[numTotal + i] = i; + widths->data.S32[numTotal + i] = 0; + locations->data.S32[numTotal - i] = - i; + widths->data.S32[numTotal - i] = 0; + } + for (int i = numInner + 1; i <= numTotal; i++) { + locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] + + widths->data.S32[numTotal + i - 1] + 1; + widths->data.S32[numTotal + i] = binning - 1; + locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - + widths->data.S32[numTotal - i + 1] - binning; + widths->data.S32[numTotal - i] = binning - 1; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i], + locations->data.S32[i] + widths->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = locations->data.S32[numTotal + i]; // Location of pixel + int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel + + for (int j = - numTotal; j <= numTotal; j++) { + int v = locations->data.S32[numTotal + j]; // Location of pixel + int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel + + // Iterate over spatial order. This loop creates the terms for + // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. + for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { + for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + kernels->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, + u, uStop, v, vStop, xOrder, yOrder); + } + } + } + } + + kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; + assert(kernels->u->data.S32[kernels->subIndex] == 0 && + kernels->v->data.S32[kernels->subIndex] == 0 && + kernels->uStop->data.S32[kernels->subIndex] == 0 && + kernels->vStop->data.S32[kernels->subIndex] == 0 && + kernels->xOrder->data.S32[kernels->subIndex] == 0 && + kernels->yOrder->data.S32[kernels->subIndex] == 0); + + psFree(locations); + psFree(widths); + + return kernels; +} + + +/// Generate FRIES kernels +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner ///< Inner radius to preserve unbinned + ) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + + int fibNum = 0; // Number of Fibonacci values + int fibLast = 1, fibTotal = 2; // Fibonacci sequence + while (fibTotal < size - inner) { + int temp = fibTotal; + fibTotal += fibLast; + fibLast = temp; + fibNum++; + } + + int numInner = inner; // Number of pixels in the inner region + int numOuter = fibNum; // Number of summed pixels in the outer region + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1) * + (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions + + psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, + size, spatialOrder); // The kernels + kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + + psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + start->data.S32[numTotal] = 0; + stop->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + start->data.S32[numTotal + i] = i; + stop->data.S32[numTotal + i] = i; + start->data.S32[numTotal - i] = -i; + stop->data.S32[numTotal - i] = -i; + } + for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal; + i++, fib = (temp = fib) + fibLast, fibLast = temp) { + start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1; + stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size); + start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i]; + stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i]; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = start->data.S32[numTotal + i]; // Location of pixel + int uStop = stop->data.S32[numTotal + i]; // Width of pixel + for (int j = - numTotal; j <= numTotal; j++) { + int v = start->data.S32[numTotal + j]; // Location of pixel + int vStop = stop->data.S32[numTotal + j]; // Width of pixel + + // Iterate over spatial order. This loop creates the terms for + // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. + for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { + for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + kernels->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, + u, uStop, v, vStop, xOrder, yOrder); + } + } + } + } + + kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; + assert(kernels->u->data.S32[kernels->subIndex] == 0 && + kernels->v->data.S32[kernels->subIndex] == 0 && + kernels->uStop->data.S32[kernels->subIndex] == 0 && + kernels->vStop->data.S32[kernels->subIndex] == 0 && + kernels->xOrder->data.S32[kernels->subIndex] == 0 && + kernels->yOrder->data.S32[kernels->subIndex] == 0); + + psFree(start); + psFree(stop); + + return kernels; +} + + pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, - const psVector *sigmas, const psVector *orders) + const psVector *sigmas, const psVector *orders, int inner, + int binning) { switch (type) { @@ -195,4 +401,8 @@ case PM_SUBTRACTION_KERNEL_ISIS: return pmSubtractionKernelsISIS(size, spatialOrder, sigmas, orders); + case PM_SUBTRACTION_KERNEL_SPAM: + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning); + case PM_SUBTRACTION_KERNEL_FRIES: + return pmSubtractionKernelsFRIES(size, spatialOrder, inner); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14019) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14106) @@ -6,13 +6,16 @@ /// Type of subtraction kernel typedef enum { - PM_SUBTRACTION_KERNEL_POIS, ///< POIS kernel --- delta functions - PM_SUBTRACTION_KERNEL_ISIS ///< ISIS kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions + PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions + PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction } pmSubtractionKernelsType; /// Kernels specification typedef struct { - pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels + pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) - psVector *sigma; ///< Width of Gaussian (for ISIS) + psVector *sigma; ///< Gaussian widths (ISIS only) + psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) psVector *xOrder, *yOrder; ///< Spatial Polynomial order (for all) int subIndex; ///< Index of kernel to be subtracted (to maintain flux conservation) @@ -44,4 +47,17 @@ ); +/// Generate SPAM kernels +pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner, ///< Inner radius to preserve unbinned + int binning ///< Kernel binning factor + ); + +/// Generate FRIES kernels +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner ///< Inner radius to preserve unbinned + ); + /// Generate a kernel of a specified type pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type @@ -49,6 +65,8 @@ int spatialOrder, ///< Order of spatial variations const psVector *sigmas, ///< Gaussian widths - const psVector *orders ///< Polynomial order of gaussians - ); + const psVector *orders, ///< Polynomial order of gaussians + int inner, ///< Inner radius to preserve unbinned + int binning ///< Kernel binning factor + ); #endif Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14019) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14106) @@ -6,5 +6,4 @@ #include -#include "pmSubtraction.h" #include "pmSubtractionStamps.h" @@ -40,23 +39,28 @@ } -psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, - float threshold, float spacing) + +psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *mask, + psMaskType maskVal, psMaskType used, float threshold, + float spacing, int border) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - if (subMask) { - PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); - PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); + if (mask) { + PS_ASSERT_IMAGE_NON_NULL(mask, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL); + } + if (used == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "The mask value for used stamps cannot be zero."); + return NULL; } PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); + PS_ASSERT_INT_NONNEGATIVE(border, NULL); + PS_ASSERT_INT_LARGER_THAN(image->numCols, 2 * border, NULL); + PS_ASSERT_INT_LARGER_THAN(image->numRows, 2 * border, NULL); - // Size of image - int numRows = image->numRows; - int numCols = image->numCols; - - // Number of stamps - int xNumStamps = numCols / spacing + 1; - int yNumStamps = numRows / spacing + 1; + maskVal |= used; // Make sure we don't get stamps we've already used + int xNumStamps = (image->numCols - 2 * border) / spacing + 1; // Number of stamps in x + int yNumStamps = (image->numRows - 2 * border) / spacing + 1; // Number of stamps in y if (stamps) { @@ -74,4 +78,7 @@ } } + // Footprint of image + int numRows = image->numRows; + int numCols = image->numCols; for (int j = 0, index = 0; j < yNumStamps; j++) { @@ -89,20 +96,39 @@ // Bounds of region to search for stamp - int yMin = j * numRows / (yNumStamps + 1); - int yMax = (j + 1) * numRows / (yNumStamps + 1) - 1; - int xMin = i * numCols / (xNumStamps + 1); - int xMax = (i + 1) * numCols / (xNumStamps + 1) - 1; - assert(yMax < image->numRows && xMax < image->numCols && - yMin >= 0 && xMin >= 0); + int yMin = border + j * (numRows - border) / (yNumStamps + 1); + int yMax = border + (j + 1) * (numRows - border) / (yNumStamps + 1) - 1; + int xMin = border + i * (numCols - border) / (xNumStamps + 1); + int xMax = border + (i + 1) * (numCols - border) / (xNumStamps + 1) - 1; + assert(yMax < image->numRows - border && xMax < image->numCols - border && + yMin >= border && xMin >= border); for (int y = yMin; y <= yMax ; y++) { for (int x = xMin; x <= xMax ; x++) { - if ((!subMask || !(subMask->data.PS_TYPE_MASK_DATA[y][x] & - (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT | - PM_SUBTRACTION_MASK_REJ))) && - image->data.F32[y][x] > fluxBest) { - fluxBest = image->data.F32[y][x]; - xBest = x; - yBest = y; + if (image->data.F32[y][x] > fluxBest) { + bool ok = true; + if (mask) { + // Check kernel footprint for bad pixels + if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) { + ok = false; + } else { + for (int v = -border; v <= border && ok; v++) { + for (int u = -border; u <= border && ok; u++) { + if (mask->data.PS_TYPE_MASK_DATA[y + v][x + u] & maskVal) { + ok = false; +#if 0 + // Mark it so we don't have to look so hard next time + mask->data.PS_TYPE_MASK_DATA[y][x] |= maskVal; +#endif + } + } + } + } + } + + if (ok) { + fluxBest = image->data.F32[y][x]; + xBest = x; + yBest = y; + } } } Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14019) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14106) @@ -3,4 +3,5 @@ #include + /// Status of stamp @@ -25,6 +26,9 @@ const psImage *image, ///< Image for which to find stamps const psImage *mask, ///< Mask + psMaskType maskVal, ///< Value for mask + psMaskType bad, ///< Mask value for bad stamps float threshold, ///< Threshold for stamps in the image - float spacing ///< Rough spacing for stamps + float spacing, ///< Rough spacing for stamps + int border ///< Size of border around image; no stamps in border );