Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 13736) +++ 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; + } }