Index: /trunk/psModules/src/imcombine/Makefile.am =================================================================== --- /trunk/psModules/src/imcombine/Makefile.am (revision 13339) +++ /trunk/psModules/src/imcombine/Makefile.am (revision 13340) @@ -4,13 +4,17 @@ libpsmodulesimcombine_la_LDFLAGS = -release $(PACKAGE_VERSION) -libpsmodulesimcombine_la_SOURCES = \ - pmReadoutCombine.c \ - pmImageCombine.c -# pmImageSubtract.c +libpsmodulesimcombine_la_SOURCES = \ + pmReadoutCombine.c \ + pmImageCombine.c \ + pmSubtraction.c \ + pmSubtractionKernels.c \ + pmSubtractionStamps.c pkginclude_HEADERS = \ - pmReadoutCombine.h \ - pmImageCombine.h -# pmImageSubtract.h + pmReadoutCombine.h \ + pmImageCombine.h \ + pmSubtraction.h \ + pmSubtractionKernels.h \ + pmSubtractionStamps.h CLEANFILES = *~ Index: unk/psModules/src/imcombine/pmImageSubtract.c =================================================================== --- /trunk/psModules/src/imcombine/pmImageSubtract.c (revision 13339) +++ (revision ) @@ -1,1393 +1,0 @@ -/** @file ImageSubtract.c - * - * This file will contain code which creates a set of kernel basis - * functions, solves for their solution, and applies them to an image. - * - * @author Paul Price, IfA (original prototype) - * @author GLG, MHPCC - * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-24 22:55:05 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - * XXX: sync with IfA on this: - * The (x, y) (row, col) issue is becoming a problem. In this file, and I - * think, the rest of psLib and psModules, the following conventions are used: - * - * 1) x will correspond to the column, and y will correspond to the row. - * 2) When used in function prototypes, the column (and hence x) appears - * first. - * 3) When used to index 2-D arrays, obviously, the row (and hence, y) - * appears first. (2 and 3 are the source of confusion). - * 4) When (u, v) are used in certain structures. - * u corresponds to x - * v corresponds to y - * 5) When element (a, b) is casually referred to in comments, or - * documentation it is unclear where a is the row, or the column. - * 6) A convention on loop index variables (i, j) would be convenient. - * Currently, sometimes i corresponds to the column (x), - * usually it corresponds to the row (y). - * - * XXX: The following variables are used an interpreted this way: - * kernelSize: Means that the actual kernel is a square (1 + 2 * kernelSize) per side. - * border: When accessing an image, a swath of pixels this wide is ignored. - * footprint: When accessing a stample, a square of pixels, footprint pixels per side, - * are looked at. We must ensure that (footprint+kernelSize) pixels exist - * around the center. - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include "pmImageSubtract.h" - -/******************************************************************************* - Private alloc/free functions. -XXX: It's not clear if the SubtractionKernels alloc/free functions are feasable. - ******************************************************************************/ -void p_pmStampFree(pmStamp *stamp) -{ - psFree(stamp->matrix); - psFree(stamp->vector); -} - -pmStamp *p_pmStampAlloc(pmStampStatus status) -{ - pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp)); - stamp->x = 0; - stamp->p_xSize = 0; - stamp->y = 0; - stamp->p_ySize = 0; - stamp->matrix = NULL; - stamp->vector = NULL; - stamp->status = status; - - psMemSetDeallocator(stamp, (psFreeFunc)p_pmStampFree); - - return(stamp); -} - -void p_pmSubtractionKernelsFree(psSubtractionKernels *kernels) -{ - psFree(kernels->u); - psFree(kernels->v); - psFree(kernels->sigma); - psFree(kernels->xOrder); - psFree(kernels->yOrder); - psFree(kernels->preCalc); - - psFree(kernels); -} - -psSubtractionKernels *p_pmSubtractionKernelsAlloc(int numBasisFunctions, - pmSubtractionKernelsType type) -{ - psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); - - tmp->type = type; - psMemSetDeallocator(tmp, (psFreeFunc) p_pmSubtractionKernelsFree); - return(tmp); -} - -/******************************************************************************* -psSubtractionKernels struct. - ******************************************************************************/ -psSubtractionKernels *pmSubtractionKernelsAllocPOIS(int size, - int spatialOrder) -{ - psTrace("psModules.imcombine", 3, - "Calling pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder); - PS_ASSERT_INT_POSITIVE(size, NULL); - PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); - // - // Calculate the number of basis functions (nBF) - // - psS32 xKernelHalfSize = size; - psS32 yKernelHalfSize = size; - psS32 nBF = (2 * xKernelHalfSize + 1) * - (2 * yKernelHalfSize + 1) * - (spatialOrder + 1) * - (spatialOrder + 2) / 2; - - // - // Generate the new psSubtractionKernels data structure: - // - psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); - tmp->type = PM_SUBTRACTION_KERNEL_POIS; - tmp->u = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->v = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->sigma = NULL; - tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->subIndex = 0; - tmp->preCalc = NULL; - tmp->size = size; - tmp->spatialOrder = spatialOrder; - - // - // This corresponds to the Kernel Basis Function with (u, v) = (0, 0) - // - // Put the (u,v) = (0,0) component right at the start of the list - // for convenience. - // - psS32 ptr = 0; - for (psS32 order = 0; order <= spatialOrder; order++) { - for (psS32 xOrder = 0; xOrder <= order; xOrder++) { - psS32 yOrder = order - xOrder; - tmp->u->data.F32[ptr] = 0; - tmp->v->data.F32[ptr] = 0; - tmp->xOrder->data.F32[ptr] = xOrder; - tmp->yOrder->data.F32[ptr] = yOrder; - ptr++; - } - } - - // - // Iterate over (u,v). Generate a set of kernels for each (u, v). - // - for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { - for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { - // Already did (u,v) = (0,0): it's at the start, so skip it now. - if ((u != 0) || (v != 0)) { - // - // Iterate over spatial order. This loop creates the terms for - // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. - // - for (psS32 order = 0; order <= spatialOrder; order++) { - for (psS32 xOrder = 0; xOrder <= order; xOrder++) { - psS32 yOrder = order - xOrder; - tmp->u->data.F32[ptr] = u; - tmp->v->data.F32[ptr] = v; - tmp->xOrder->data.F32[ptr] = xOrder; - tmp->yOrder->data.F32[ptr] = yOrder; - ptr++; - } - } - } - } - } - - psTrace("psModules.imcombine", 3, - "Exiting pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder); - return(tmp); -} - -/******************************************************************************* -XXX: Get the types correct (u, v, xOrder, yOrder). - -XXX: Code review this. - ******************************************************************************/ -psSubtractionKernels *pmSubtractionKernelsAllocISIS(const psVector *sigmas, - const psVector *orders, - int size, - int spatialOrder) -{ - PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); - PS_ASSERT_VECTOR_NON_NULL(orders, NULL); - psTrace("psModules.imcombine", 3, - "Calling pmSubtractionKernelsAllocISIS(%ld, %ld, %d, %d)\n", - sigmas->n, orders->n, size, spatialOrder); - PS_ASSERT_INT_POSITIVE(size, NULL); - PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); - PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); - // - // Calculate the number of basis functions (nBF). - // - psS32 numSigmas = sigmas->n; - - // XXX: Get rid of the sigma loop? We merely multiple nBF by numSigmas? - // XXX: Verify that all this is correct. - psS32 nBF = 0; - for (psS32 s = 0 ; s < numSigmas ; s++) { - for (psS32 o = 0 ; o < orders->n ; o++) { - nBF+=((orders->data.S32[o] + 1) * (orders->data.S32[o] + 2) / 2); - } - } - nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2); - - // - // Generate the new psSubtractionKernels data structure: - // - psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); - tmp->type = PM_SUBTRACTION_KERNEL_ISIS; - tmp->u = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->v = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->sigma = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->subIndex = 0; - tmp->size = size; - tmp->spatialOrder = spatialOrder; - tmp->preCalc = psArrayAlloc(nBF); - - // - // We loop through all combinations of sigmas and polynomial orders - // creating the kernel basis functions. - // - psS32 ptr = 0; - for (psS32 sigPtr = 0 ; sigPtr < numSigmas ; sigPtr++) { - tmp->sigma->data.F32[sigPtr] = sigmas->data.F32[sigPtr]; - // - // (xOrderP, yOrderP) are the order of the polynomial that modify the - // gaussian in the kernel. They correspond to (j, k) in equation (5) - // from the psModules SDRS. - // - for (psS32 o = 0 ; o < orders->n ; o++) { - for (psS32 orderP = 0 ; orderP <= orders->data.S32[o] ; orderP++) { - for (psS32 xOrderP = 0 ; xOrderP <= orderP ; xOrderP++) { - psS32 yOrderP = orderP - xOrderP; - - psImage *currPreCalc = psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); - psImageInit(currPreCalc, 0.0); - psBool setPreCalc = true; - // - // We loop through all spatial orders. Since they have no effect on - // the preCalc images, we only calculate them once, and store pointers - // in tmp->preCalc->data[ptr] for other spatial orders. - // - for (psS32 order = 0; order <= spatialOrder; order++) { - for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { - PS_ASSERT_INT_LESS_THAN(ptr, nBF, NULL); - - psS32 orderYTerm = order - orderXTerm; - - tmp->u->data.F32[ptr] = xOrderP; - tmp->v->data.F32[ptr] = yOrderP; - tmp->xOrder->data.F32[ptr] = orderXTerm; - tmp->yOrder->data.F32[ptr] = orderYTerm; - tmp->sigma->data.F32[ptr] = sigmas->data.F32[sigPtr]; - tmp->preCalc->data[ptr] = (psPtr *) currPreCalc; - - // - // We calculate the preCalc image only the first time through - // this loop. Otherwise, we increment the memory reference - // counter. - // - if (setPreCalc == true) { - for (psS32 v = -size; v <= size; v++) { - for (psS32 u = -size; u <= size; u++) { - // Scale the (u,v) coordinates in kernel space to [-1.0:1.0]. - psF32 uScaled = ((psF32) u) / ((psF32) size); - psF32 vScaled = ((psF32) v) / ((psF32) size); - - // Compute the value of the kernel at location (u, v): - psF32 exponent = (PS_SQR(uScaled) + PS_SQR(vScaled)) / - (2.0 * PS_SQR(sigmas->data.F32[sigPtr])); - currPreCalc->data.F32[v+size][u+size] = - exp(-exponent) * - pow(uScaled, orderXTerm) * - pow(vScaled, orderYTerm); - } - } - setPreCalc = false; - } else { - psMemIncrRefCounter(currPreCalc); - } - ptr++; - } - } - } - } - } - } - - psTrace("psModules.imcombine", 3, - "Exiting pmSubtractionKernelsAllocISIS(%ld, %ld, %d, %d)\n", - sigmas->n, orders->n, size, spatialOrder); - return(tmp); -} - -/******************************************************************************* -pmSubtractionFindStamps(stamps, image, mask, maskVal, threshold, xNum, yNum, border) - -XXX: The SDRS and the proptotype code differ significantly. - Prototype: When a maximum pixel is found within a stamp, an area of size - 2*footprint is searched around that pixel looking for masked pixels. - SDRS: none of that is required. - -XXX: Do we need to care about case where yNum/xNum does not evenly divide the -nnumber of rows/columns in the image? - ******************************************************************************/ -psArray *pmSubtractionFindStamps(psArray *stamps, ///< Output stamps, or NULL - const psImage *image, ///< Image for which to find stamps - const psImage *mask, ///< Mask - psU32 maskVal, ///< Value for mask - psF32 threshold, ///< Threshold for stamps in the image - psS32 xNum, ///< Number of stamps in x - psS32 yNum, ///< Number of stamps in y - psS32 border ///< Border around image to ignore (should be size of kernel or larger) - ) -{ - psTrace("psModules.imcombine", 3, - "Calling pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n", - maskVal, threshold, xNum, yNum, border); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - if (mask != NULL) { - PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_INT_POSITIVE(xNum, NULL); - PS_ASSERT_INT_POSITIVE(yNum, NULL); - PS_ASSERT_INT_POSITIVE(border, NULL); - PS_ASSERT_INT_LARGER_THAN(image->numCols, xNum, NULL); - PS_ASSERT_INT_LARGER_THAN(image->numRows, yNum, NULL); - PS_ASSERT_INT_LARGER_THAN(image->numCols, (2 * border), NULL); - PS_ASSERT_INT_LARGER_THAN(image->numRows, (2 * border), NULL); - - if (stamps != NULL) { - PS_ASSERT_INT_EQUAL(stamps->n, (xNum * yNum), NULL); - // - // Ensure that a pmStamp struct exists at each psArray location. - // - for (psS32 s = 0 ; s < (xNum * yNum) ; s++) { - if (NULL == stamps->data[s]) { - stamps->data[s] = (psPtr *) p_pmStampAlloc(PM_STAMP_REJECTED); - } - } - } else { - stamps = (psArray *) psArrayAlloc(xNum * yNum); - for (psS32 s = 0 ; s < (xNum * yNum) ; s++) { - stamps->data[s] = (psPtr *) p_pmStampAlloc(PM_STAMP_REJECTED); - } - } - psS32 numRows = image->numRows; - psS32 numCols = image->numCols; - - // - // Iterate over the image sections - // - // XXX: Must handle cases where image size is not an even multiple of xNum or yNum - // they are currently ignored. - // - psS32 num = 0; - for (psS32 j = 0; j < yNum; j++) { - for (psS32 i = 0; i < xNum; i++) { - pmStamp *stamp = (pmStamp *) stamps->data[num]; - // - // Only find a new stamp if we need to - // - if (stamp->status == PM_STAMP_REJECTED) { - // - // Find maximum non-masked value in the image section, - // but don't include a footprint around the edge - // - psF32 max = -INFINITY; - psS32 bestx = 0; - psS32 besty = 0; - // - // The following nested loop iterates over every pixel in the mask - // associated with this (i, j). It ignores pixels within a - // border of pixels from the image edge. - // - // XXX: verify (numX, numY), then get rid of it. - // - psS32 numX = xNum; - psS32 numY = yNum; - psS32 yMin = border + j * (numCols - 2.0 * border) / numY; - psS32 yMax = (border + (j + 1) * (numCols - 2.0 * border) / numY) - 1; - psS32 xMin = border + i * (numRows - 2.0 * border) / numX; - psS32 xMax = (border + (i + 1) * (numRows - 2.0 * border) / numX) - 1; - - if ((yMax >= image->numRows) || - (xMax >= image->numCols) || - (yMin < 0) || - (xMin < 0)) { - // XXX: We skip this stamp since its borders extends beyond the image. - // XXX: This is here mainly as a safeguard. We need to redefine the above - // min/max pixels calculation to ensure that all stamps are legitimate. - - stamp->x = -1; - stamp->y = -1; - stamp->status = PM_STAMP_NONE; - } else { - stamp->p_xSize = 1 + (xMax - xMin); - stamp->p_ySize = 1 + (yMax - yMin); - stamp->p_xMin = xMin; - stamp->p_xMax = xMax; - stamp->p_yMin = yMin; - stamp->p_yMax = yMax; - - for (psS32 y = yMin; y <= yMax ; y++) { - for (psS32 x = xMin; x <= xMax ; x++) { - // Determine if this pixel is larger than the max, and unmasked. - if (image->data.F32[y][x] > max) { - if ((mask == NULL) || !((mask->data.U8[y][x]) & maskVal)) { - max = image->data.F32[y][x]; - bestx = x; - besty = y; - } - } - } - } - - // - // If the max pixel is larger than the threshold, we keep this stamp. - // Otherwise, mark the stamp as PM_STAMP_NONE - // - if (image->data.F32[besty][bestx] >= threshold) { - stamp->x = bestx; - stamp->y = besty; - stamp->status = PM_STAMP_RECALC; - } else { - stamp->x = bestx; - stamp->y = besty; - stamp->status = PM_STAMP_NONE; - } - } - } - num++; - } - } - psTrace("psModules.imcombine", 3, - "Exiting pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n", - maskVal, threshold, xNum, yNum, border); - return(stamps); -} - -/******************************************************************************* -GenSpatialOrder(spatialOrder, x, y): generates and returns a psImage in which -the [i][j] location is calculated as (x^i * y^j). - -XXX: Modify loop so that terms higher than spatialOrder are not computed. - -XXX: Modify this so that [i][j] location is calculated as (x^j * y^i)? - ******************************************************************************/ -static psImage *GenSpatialOrder(psS32 spatialOrder, - psF32 x, - psF32 y) -{ - psTrace("psModules.imcombine", 4, - "Calling GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y); - - psImage *polyValues = psImageAlloc(spatialOrder+1, spatialOrder+1, PS_TYPE_F64); - - psF64 xSum = 1.0; - psF64 ySum = 1.0; - for (psS32 i = 0; i < spatialOrder + 1; i++) { - ySum = xSum; - for (psS32 j = 0; j < spatialOrder + 1; j++) { - polyValues->data.F64[i][j] = ySum; - ySum*= y; - } - xSum*= x; - } - - psTrace("psModules.imcombine", 4, - "Exiting GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y); - - return(polyValues); -} - - -/******************************************************************************* -IsisKernelConvolve(input, kernels, kernelID, col, row): This routine -convolves a single kernel basis function with a pixel in an image. - ******************************************************************************/ -static psF32 IsisKernelConvolve(const psImage *input, - const psSubtractionKernels *kernels, - psS32 kernelID, - psS32 col, - psS32 row) -{ - - - psTrace("psModules.imcombine", 4, - "Calling IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row); - psS32 spatialOrder = kernels->spatialOrder; - psS32 kernelSize = kernels->size; - psS32 xOrder = (psS32) kernels->xOrder->data.F32[kernelID]; - psS32 yOrder = (psS32) kernels->yOrder->data.F32[kernelID]; - psF32 numColsHalf = 0.5 * (psF32) input->numCols; - psF32 numRowsHalf = 0.5 * (psF32) input->numRows; - psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position - psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position - - psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); - - psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; - - psImage *preCalc = (psImage *) kernels->preCalc->data[kernelID]; - - // XXX: Are the following asserts really necessary? - PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(row-kernelSize, 0, NAN); - PS_ASSERT_INT_LESS_THAN(row+kernelSize, input->numRows, NAN); - PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(col-kernelSize, 0, NAN); - PS_ASSERT_INT_LESS_THAN(col+kernelSize, input->numCols, NAN); - psF32 conv = 0.0; - for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { - for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { - conv += input->data.F32[yy+row][xx+col] * - preCalc->data.F32[yy+kernelSize][xx+kernelSize] * - polyVal; - } - } - psFree(polyValues); - - psTrace("psModules.imcombine", 4, - "Exiting IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row); - return(conv); -} - -/******************************************************************************* -ConvolvePixelPois(input, mask, badStampMaskVal, solution, kernels, col, row): - -This routine takes a single pixel in the psImage input and convolves it with -the set of kernel basis functions and their appropriate weights in solution. -It returns the value of the convolved pixel. - -XXX: Static structure for polyValues? - ******************************************************************************/ -static psF32 ConvolvePixelPois(const psImage *input, - const psImage *mask, - psU32 badStampMaskVal, - const psVector *solution, - const psSubtractionKernels *kernels, - psS32 col, - psS32 row) -{ - psTrace("psModules.imcombine", 4, - "Calling ConvolvePixelPois(%d, %d)\n", col, row); - psS32 nBF = kernels->u->n; - psF32 numColsHalf = 0.5 * (psF32) input->numCols; - psF32 numRowsHalf = 0.5 * (psF32) input->numRows; - psF32 background = solution->data.F64[solution->n-1]; - psS32 spatialOrder = kernels->spatialOrder; - psF32 conv = background; // Initial convolved value - - if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { - psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position - psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position - psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); - - // Iterate over the kernel basis functions - for (psS32 k = 0; k < nBF; k++) { - psS32 u = (psS32) kernels->u->data.F32[k]; - psS32 v = (psS32) kernels->v->data.F32[k]; - - // XXX: What's the story with this? - #if 0 - - psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; - psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; - psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; - #else - - psF32 polyVal = 1.0; - #endif - - // XXX: Why this? - if (k == 0) { - conv += solution->data.F64[k] * input->data.F32[row - v][col - u] * polyVal; - } else { - conv += solution->data.F64[k] * - (input->data.F32[row - v][col - u] * polyVal - input->data.F32[row][col]); - } - } - psFree(polyValues); - } - - psTrace("psModules.imcombine", 4, - "Exiting ConvolvePixelPois(%d, %d)\n", col, row); - return(conv); -} - - - -/******************************************************************************* -ConvolvePixelIsis(input, mask, badStampMaskVal, solution, kernels, col, row): - -This routine takes a single pixel in the psImage input and convolves it with -the set of kernel basis functions and their appropriate weights in solution. -It returns the value of the convolved pixel. - -XXX: Static structure for polyValues? - ******************************************************************************/ -static psF32 ConvolvePixelIsis(const psImage *input, - const psImage *mask, - psU32 badStampMaskVal, - const psVector *solution, - const psSubtractionKernels *kernels, - psS32 col, - psS32 row) -{ - psTrace("psModules.imcombine", 4, - "Calling ConvolvePixelIsis(%d, %d)\n", col, row); - psF32 background = solution->data.F64[solution->n-1]; - psF32 conv = background; // Initial convolved value - - if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { - // Iterate over the kernel basis functions - for (psS32 k = 0; k < kernels->u->n; k++) { - conv += IsisKernelConvolve(input, kernels, k, col, row); - } - } - - psTrace("psModules.imcombine", 4, - "Exiting ConvolvePixelIsis(%d, %d)\n", col, row); - return(conv); -} - -/******************************************************************************* -ConvolveImage(input, mask, badStampMaskVal, solution, kernels): convolves an -arbitrary image with either an ISIS or POIS set of kernel basis functions. - ******************************************************************************/ -static psImage *ConvolveImage(const psImage *input, - const psImage *mask, - psU32 badStampMaskVal, - const psVector *solution, - const psSubtractionKernels *kernels) -{ - psTrace("psModules.imcombine", 4, "Calling ConvolveImage()\n"); - PS_ASSERT_IMAGE_NON_NULL(input, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(input, NULL); - PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, NULL); - if (mask != NULL) { - PS_ASSERT_IMAGES_SIZE_EQUAL(input, mask, NULL); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(solution, NULL); - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - if (kernels->preCalc != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); - } else { - if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); - return(NULL); - } - } - psS32 nBF = kernels->u->n; - PS_ASSERT_VECTOR_SIZE(solution, (long)nBF+1, NULL); - - psS32 numCols = input->numCols; - psS32 numRows = input->numRows; - psS32 kernelSize = kernels->size; - - psImage *convolved = psImageAlloc(numCols, numRows, PS_TYPE_F32); - - for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { - for (psS32 x = kernelSize; x < numCols - kernelSize; x++) { - if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { - convolved->data.F32[y][x] = ConvolvePixelPois(input, mask, badStampMaskVal, - solution, kernels, x, y); - } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - convolved->data.F32[y][x] = ConvolvePixelIsis(input, mask, badStampMaskVal, - solution, kernels, x, y); - } else { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel type. Returning NULL\n"); - return(NULL); - } - } - } - - // - // Pad the rest of the convolved image with 0.0 - // - for (psS32 y = kernelSize; y < numRows - kernelSize; y++) { - for (psS32 x = 0; x < kernelSize; x++) { - convolved->data.F32[y][x] = 0.0; - } - for (psS32 x = numCols - kernelSize; x < numCols; x++) { - convolved->data.F32[y][x] = 0.0; - } - } - for (psS32 y = 0; y < kernelSize; y++) { - for (psS32 x = 0; x < numCols; x++) { - convolved->data.F32[y][x] = 0.0; - } - } - for (psS32 y = numRows - kernelSize; y < numRows; y++) { - for (psS32 x = 0; x < numCols; x++) { - convolved->data.F32[y][x] = 0.0; - } - } - - psTrace("psModules.imcombine", 4, "Exiting ConvolveImage()\n"); - return convolved; -} - - - - - -/******************************************************************************* -XXX: We should assert that the (footprint, kernelSize, imageSize) stuff -ensures that all data is accessed in bounds? - ******************************************************************************/ -bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, - const psImage *reference, ///< Reference image - const psImage *input, ///< Input image - const psSubtractionKernels *kernels, ///< The kernel basis functions - psS32 footprint ///< Half-size of region over which to calculate equation - ) -{ - psTrace("psModules.imcombine", 3, - "Calling pmSubtractionCalculateEquation()\n"); - PS_ASSERT_PTR_NON_NULL(stamps, false); - PS_ASSERT_IMAGE_NON_NULL(reference, false); - PS_ASSERT_IMAGE_NON_EMPTY(reference, false); - PS_ASSERT_IMAGE_NON_NULL(input, false); - PS_ASSERT_IMAGE_NON_EMPTY(input, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference, input, false); - PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); - if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); - } - psS32 kernelSize = kernels->size; - PS_ASSERT_INT_NONNEGATIVE(footprint, false); - // - // For each legitimate stamp, ensure that the footprint is small enough to perform - // the full calculation. - // - // XXX: Verify with IfA that this is a reasonable action. - // - for (psS32 s = 0; s < stamps->n; s++) { - pmStamp *stamp = (pmStamp *) stamps->data[s]; - if (stamp->status == PM_STAMP_RECALC) { - // XXX: trace message - // printf("stamp %d (x, y) is (%d, %d). Footprint is %d. kernelSize is %d.\n", s, stamp->x, stamp->y, footprint, kernelSize); - if (((stamp->y - (footprint + kernelSize)) < 0) || - ((stamp->x - (footprint + kernelSize)) < 0) || - ((stamp->y + footprint + kernelSize) >= input->numRows) || - ((stamp->x + footprint + kernelSize) >= input->numCols)) { - stamp->status = PM_STAMP_NONE; - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: stamp %d will be ignored. It exceeds image size: access columns (%d to %d) and rows (%d to %d)\n", - s, - stamp->x - (footprint + kernelSize), - (stamp->x + footprint + kernelSize) - 1, - stamp->y - (footprint + kernelSize), - (stamp->y + footprint + kernelSize) - 1); - } - } - } - - psS32 numHalfRows = reference->numRows; - psS32 numHalfCols = reference->numCols; - psS32 spatialOrder = kernels->spatialOrder; - - // - // The numSolveParams incorporates the additional parameter for the - // background value, which we must solve for. - // - psS32 numKernels = kernels->u->n; - int numSolveParams = numKernels + 1; - int bgIndex = numKernels; // Index in matrix for the background - - // - // We iterate over each stamp, allocate the matrix and vectors if - // necessary, and then calculate those matrix/vectors. - // - for (psS32 s = 0; s < stamps->n; s++) { - pmStamp *stamp = (pmStamp *) stamps->data[s]; - psTrace("psModules.imcombine", 6, "subCalcEqn(): stamp %d\n", s); - if (stamp->status == PM_STAMP_RECALC) { - psTrace("psModules.imcombine", 6, "subCalcEqn(): stamp %d: status is PM_STAMP_RECALC.\n", s); - psImage *stampMatrix = stamp->matrix; - psVector *stampVector = stamp->vector; - - if (stampMatrix == NULL) { - stampMatrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); - stamp->matrix = stampMatrix; - } else { - PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, false); - PS_ASSERT_IMAGE_SIZE(stampMatrix, numSolveParams, numSolveParams, false); - } - psImageInit(stampMatrix, 0.0); - - if (stampVector == NULL) { - stampVector = psVectorAlloc(numSolveParams, PS_TYPE_F64); - stamp->vector = stampVector; - } else { - PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, false); - PS_ASSERT_VECTOR_SIZE(stampVector, (long)numSolveParams, false); - } - psVectorInit(stampVector, 0.0); - psTrace("psModules.imcombine", 6, "subCalcEqn(): stamp %d: allocate matrix and vector.\n", s); - - // - // Evaluate the spatial-order polynomial. The [i][j]-th element of - // the psImage polyValues will hold (x^i * y^j) for the stamp. The - // (x, y) value are scaled to [-1:1] - // - psImage *polyValues = GenSpatialOrder(spatialOrder, - ((psF64) (stamp->x - numHalfCols)) / ((psF64) numHalfCols), - ((psF64) (stamp->y - numHalfRows)) / ((psF64) numHalfRows)); - - psTrace("psModules.imcombine", 6, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s); - - if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { - // - // Iterate over all pixels surrounding this stamp. - // - for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { - for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { - psTrace("psModules.imcombine", 6, "subCalcEqn(): pixel (%d, %d).\n", y, x); - - // The inverse of the noise, squared. - psF32 invNoise2 = 1.0/reference->data.F32[y][x]; - - // - // Iterate over the first convolution */ - // - for (psS32 k1 = 0; k1 < numKernels; k1++) { - psS32 u1 = kernels->u->data.F32[k1]; // Offset in x - psS32 v1 = kernels->v->data.F32[k1]; // Offset in y - psS32 i1 = kernels->xOrder->data.F32[k1]; // Polynomial order in x - psS32 j1 = kernels->yOrder->data.F32[k1]; // Polynomial order in y - - // - // First convolution. This will set the value for the stampVector. - // - // XXX: verify the [y-v2][x-u2] subscript. This generated errors in - // testing, depending on kernel size and footprint. - // - psF32 conv1 = polyValues->data.F64[j1][i1] * reference->data.F32[y - v1][x - u1]; - - // - // Assuming that the first kernel component is 0 order in x and y, and 0 offset - // XXX: I don't understand this: - // - if (k1 != 0) { - conv1 -= reference->data.F32[y][x]; - } - - // - // Iterate over the second convolution - // - for (psS32 k2 = k1; k2 < numKernels; k2++) { - psS32 u2 = (psS32) kernels->u->data.F32[k2]; // Offset in x - psS32 v2 = (psS32) kernels->v->data.F32[k2]; // Offset in y - psS32 i2 = (psS32) kernels->xOrder->data.F32[k2]; // Polynomial order in x - psS32 j2 = (psS32) kernels->yOrder->data.F32[k2]; // Polynomial order in y - // - // XXX: verify the [y-v2][x-u2] subscript. This generated errors in - // testing, depending on kernel size and footprint. - // - // Second convolution - // - psF32 conv2 = polyValues->data.F64[j2][i2] * - reference->data.F32[y-v2][x-u2]; - // - // Assuming that the first kernel component is 0 order in x and y, and 0 offset - // XXX: I don't understand this: - if (k2 != 0) { - // - conv2 -= reference->data.F32[y][x]; - } - - // Add into the matrix element - stampMatrix->data.F64[k1][k2] += conv1 * conv2 * invNoise2; - - } // Iteration on second convolution - - // Add into the vector element - stampVector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2; - - /* Background term */ - stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2; - - } // Iteration on first convolution - - // - // Background only terms. - // XXX: understand this. - // - stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; - stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; - } - } - } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - // - // Iterate over all pixels surrounding this stamp. - // - // XXX: Why isn't there a polyValues term here? - // - - for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { - for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { - psTrace("psModules.imcombine", 6, "subCalcEqn(): pixel (%d, %d).\n", y, x); - psF32 invNoise2 = 1.0/reference->data.F32[y][x]; // The inverse of the noise, squared. - - for (psS32 k1 = 0; k1 < numKernels; k1++) { - psF32 conv1 = IsisKernelConvolve(reference, kernels, k1, x, y); - - for (psS32 k2 = k1; k2 < numKernels; k2++) { - //printf("(k1, k2) is (%d, %d)\n", k1, k2); - psF32 conv2 = IsisKernelConvolve(reference, kernels, k2, x, y); - stampMatrix->data.F64[k1][k2]+= conv1 * conv2 * invNoise2; - } - stampVector->data.F64[k1]+= input->data.F32[y][x] * conv1 * invNoise2; - stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2; - } - stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; - stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; - } - } - } else { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel->type.\n"); - return(false); - } - psFree(polyValues); - - // XXX: Generate psTrace() - if (0) { - for (psS32 s = 0; s < stamps->n; s++) { - pmStamp *stamp = (pmStamp *) stamps->data[s]; - if (stamp->status == PM_STAMP_RECALC) { - psVector *stampVector = stamp->vector; - printf("STAMP: stamp %d vector:\n", s); - PS_VECTOR_PRINT_F64(stampVector); - } - } - } - - // - // Fill in other side of symmetric matrix - // - // XXX: understand this. - // XXX: Why aren't they using numSolveParams instead of numKernels? - // XXX: is this POIS specific? - // - for (psS32 k1 = 0; k1 < numKernels; k1++) { - for (psS32 k2 = 0; k2 < k1; k2++) { - stampMatrix->data.F64[k1][k2] = stampMatrix->data.F64[k2][k1]; - } - stampMatrix->data.F64[bgIndex][k1] = stampMatrix->data.F64[k1][bgIndex]; - } - - // - // XXX: Why aren't they using numSolveParams instead of numKernels? - // XXX: is this POIS specific? - // - #define XXX_CONFIG_PENALTY 1.0 - for (psS32 k = 0; k < numKernels; k++) - { - psS32 u = kernels->u->data.F32[k]; // Offset in x - psS32 v = kernels->v->data.F32[k]; // Offset in y - stampMatrix->data.F64[k][k] += XXX_CONFIG_PENALTY * (psF32)(u*u + v*v); - } - stamp->status = PM_STAMP_USED; - } else { - // Stamp is ignored since it's not PM_STAMP_RECALC - } - } - psTrace("psModules.imcombine", 3, - "Exiting pmSubtractionCalculateEquation()\n"); - return(true); -} - - - - -/******************************************************************************* - ******************************************************************************/ -psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL - const psArray *stamps ///< Array of stamps - ) -{ - psTrace("psModules.imcombine", 3, - "Calling pmSubtractionSolveEquation()\n"); - PS_ASSERT_PTR_NON_NULL(stamps, NULL); - psS32 size = -1; - psS32 s = 0; - - // - // Determine the size of the stamp vectors and matrix. - // We iterate until we find the first acceptable stamp. - // - while ((size == -1) && (s < stamps->n)) { - pmStamp *stamp = (pmStamp *) stamps->data[s]; - PS_ASSERT_PTR_NON_NULL(stamp, NULL); - if (stamp->status == PM_STAMP_USED) { - size = ((pmStamp *) stamps->data[s])->vector->n; - PS_ASSERT_INT_POSITIVE(size, NULL); - } - s++; - } - if (size == -1) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: no acceptable stamps. Returning NULL\n"); - return(NULL); - } - - if (solution != NULL) { - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_SIZE(solution, (long)size, NULL); - } else { - solution = psVectorAlloc(size, PS_TYPE_F64); - } - - // - // Create the solution matrix and vector. - // - // XXX: Test these functions with size=-1. This caused seg faults during test. - // This should be done in the psImage.c and psVector.c test files. It - // should never occur here. - // - psImage *sumMatrix = psImageAlloc(size, size, PS_TYPE_F64); - psVector *sumVector = psVectorAlloc(size, PS_TYPE_F64); - psVectorInit(sumVector, 0.0); - psImageInit(sumMatrix, 0.0); - - // - // Verify that all stamps have similar sizes. - // Compute the sum matrix and vector. - // - for (psS32 s = 0; s < stamps->n; s++) { - pmStamp *stamp = (pmStamp *) stamps->data[s]; - - if (stamp->status == PM_STAMP_USED) { - PS_ASSERT_INT_EQUAL(((pmStamp *) stamps->data[s])->vector->n, size, NULL); - - psImage *stampMatrix = stamp->matrix; - psVector *stampVector = stamp->vector; - PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_SIZE(stampVector, (long)size, NULL); - PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, NULL); - PS_ASSERT_IMAGE_SIZE(stampMatrix, size, size, NULL); - - (void)psBinaryOp(sumMatrix, sumMatrix, "+", stampMatrix); - (void)psBinaryOp(sumVector, sumVector, "+", stampVector); - } - } - psVector *permutation = NULL; - // XXX: Check output from these routines. - - // XXX: psTrace() - if (0) { - PS_IMAGE_PRINT_F64(sumMatrix); - } - - psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); - if (luMatrix == NULL) { - psError(PS_ERR_UNKNOWN, true, "Failed to LU-Decompose the matrix.\n"); - psFree(sumMatrix); - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - return(NULL); - } - // XXX: psTrace() - if (0) { - PS_IMAGE_PRINT_F64(luMatrix); - } - - solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation); - // XXX: psTrace() - // XXX: should we be checking for NAN's in the solution vector? - if (0) { - PS_VECTOR_PRINT_F64(solution); - } - if (solution == NULL) { - psError(PS_ERR_UNKNOWN, true, "Failed to solve the matrix.\n"); - psFree(sumMatrix); - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - return(NULL); - } - - psFree(sumMatrix); - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - - psTrace("psModules.imcombine", 3, - "Exiting pmSubtractionSolveEquation()\n"); - return(solution); -} - - -/******************************************************************************* - ******************************************************************************/ -static psVector *CalculateDeviations(psVector *deviations, - psArray *stamps, - psS32 footprint, - const psImage *refImage, - const psImage *inImage, - const psImage *mask, - psU32 badStampMaskVal, - const psSubtractionKernels *kernels, - const psVector *solution) -{ - psTrace("psModules.imcombine", 4, - "Calling CalculateDeviations()\n"); - PS_ASSERT_PTR_NON_NULL(stamps, NULL); - if (deviations != NULL) { - PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(deviations, stamps, NULL); - } else { - deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); - // XXX: Probably not necessary. - psVectorInit(deviations, 0.0); - } - PS_ASSERT_IMAGE_NON_NULL(refImage, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(refImage, NULL); - PS_ASSERT_IMAGE_TYPE(refImage, PS_TYPE_F32, NULL); - PS_ASSERT_IMAGE_NON_NULL(inImage, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(inImage, NULL); - PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, NULL); - PS_ASSERT_IMAGE_NON_NULL(mask, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(mask, NULL); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, NULL); - PS_ASSERT_VECTOR_NON_NULL(solution, NULL); - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - if (kernels->preCalc != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); - } else { - if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); - return(NULL); - } - } - psS32 nBF = kernels->u->n; - PS_ASSERT_VECTOR_SIZE(solution, (long)nBF+1, NULL); - - psS32 kernelSize = kernels->size; - int xSize = footprint + kernelSize; - int ySize = footprint + kernelSize; - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics - psImage *subStamp = psImageAlloc(2 * xSize, 2 * ySize, PS_TYPE_F32); // Subtraction of stamp - for (psS32 s = 0; s < stamps->n; s++) { - pmStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest - psS32 x = stamp->x; // Stamp x coord - psS32 y = stamp->y; // Stamp y coord - if (stamp->status == PM_STAMP_USED) { - - psRegion myReg = psRegionSet(x - xSize, x + xSize, y - ySize, y + ySize); - psImage *refStamp = psImageSubset((psImage *) refImage, myReg); - psImage *inStamp = psImageSubset((psImage *) inImage, myReg); - psImage *maskStamp = psImageSubset((psImage *) mask, myReg); - psImage *convRefStamp = ConvolveImage(refStamp, maskStamp, badStampMaskVal, solution, kernels); - - // Calculate chi^2 - (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp); - (void)psBinaryOp(subStamp, subStamp, "/", inStamp); - (void)psBinaryOp(subStamp, subStamp, "*", subStamp); - myReg = psRegionSet(kernelSize, kernelSize + 2 * footprint, - kernelSize, kernelSize + 2 * footprint); - psImage *subStampTrim = psImageSubset((psImage *) subStamp, myReg); - psImage *maskStampTrim = psImageSubset((psImage *) maskStamp, myReg); - psImageStats(stats, subStampTrim, maskStampTrim, badStampMaskVal); - - deviations->data.F32[s] = stats->sampleMean * (psF32)footprint * (psF32)footprint * 4.0; - // XXX: Allocate and free these elsewhere. - psFree(refStamp); - psFree(inStamp); - psFree(maskStamp); - psFree(convRefStamp); - psFree(subStampTrim); - psFree(maskStampTrim); - } - } - - psFree(stats); - psFree(subStamp); - - psTrace("psModules.imcombine", 4, - "Exiting CalculateDeviations()\n"); - return deviations; -} - -/******************************************************************************* - ******************************************************************************/ -bool pmSubtractionRejectStamps(psArray *stamps, ///< Array of stamps to check for rejection - psImage *mask, ///< Mask image - psU32 badStampMaskVal, ///< Value to use in mask for bad stamp - psS32 footprint, ///< Region to mask if stamp is bad - psF32 sigmaRej, ///< Number of RMS deviations above zero at which to reject - const psImage *refImage, ///< Reference image - const psImage *inImage, ///< Input image - const psVector *solution, ///< Solution vector - const psSubtractionKernels *kernels ///< Array of kernel parameters - ) -{ - psTrace("psModules.imcombine", 3, - "Calling pmSubtractionRejectStamps()\n"); - PS_ASSERT_PTR_NON_NULL(stamps, false); - PS_ASSERT_IMAGE_NON_NULL(refImage, false); - PS_ASSERT_IMAGE_NON_EMPTY(refImage, false); - PS_ASSERT_IMAGE_TYPE(refImage, PS_TYPE_F32, false); - PS_ASSERT_IMAGE_NON_NULL(inImage, false); - PS_ASSERT_IMAGE_NON_EMPTY(inImage, false); - PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, false); - PS_ASSERT_IMAGE_NON_NULL(mask, false); - PS_ASSERT_IMAGE_NON_EMPTY(mask, false); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, false); - PS_ASSERT_VECTOR_NON_NULL(solution, false); - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false); - PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); - if (kernels->preCalc != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); - } else { - if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); - return(false); - } - } - - psS32 nBF = kernels->u->n; - PS_ASSERT_VECTOR_SIZE(solution, (long)nBF+1, false); - - psVector *deviations = CalculateDeviations(NULL, - stamps, - footprint, - refImage, - inImage, - mask, - badStampMaskVal, - kernels, - solution); - // - // Calculate the deviation from zero. - // - psF64 meanDev = 0.0; - psS32 numDev = 0; - for (psS32 i = 0; i < deviations->n; i++) { - pmStamp *stamp = stamps->data[i]; - if (stamp->status == PM_STAMP_USED) { - meanDev += PS_SQR(deviations->data.F32[i]); - numDev++; - } - } - psF32 rmsDev = sqrtf(meanDev / (psF64)(numDev - 1)); - psF32 limit = rmsDev * sigmaRej; - - for (psS32 s = 0; s < stamps->n; s++) { - pmStamp *stamp = (pmStamp *) stamps->data[s]; - if (stamp->status == PM_STAMP_USED && fabsf(deviations->data.F32[s]) > limit) { - // Mask out the stamp in the image so you don't find it again - for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { - for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { - mask->data.U8[y][x] |= badStampMaskVal; - } - } - - // Set stamp for replacement - stamp->x = 0; - stamp->y = 0; - stamp->status = PM_STAMP_REJECTED; - } - } - - psFree(deviations); - psTrace("psModules.imcombine", 3, - "Exiting pmSubtractionRejectStamps()\n"); - return(true); -} - -/******************************************************************************* - ******************************************************************************/ -psImage *pmSubtractionKernelImage(psImage *out, - const psVector *solution, - const psSubtractionKernels *kernels, - psF32 x, - psF32 y - ) -{ - psTrace("psModules.imcombine", 3, - "Calling pmSubtractionKernelImage()\n"); - PS_ASSERT_VECTOR_NON_NULL(solution, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); - if (kernels->preCalc != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); - } else { - if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); - return(NULL); - } - } - PS_ASSERT_INT_EQUAL(1+kernels->u->n, solution->n, NULL); - - psS32 nBF = kernels->u->n; - psS32 spatialOrder = kernels->spatialOrder; - psS32 kernelSize = kernels->size; - - if (out != NULL) { - if ((out->numCols < (1+2*kernelSize)) || (out->numRows < (1+2*kernelSize))) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: out image is not large enough.\n"); - return(out); - } - } else { - out = psImageAlloc(1+2*kernelSize, 1+2*kernelSize, PS_TYPE_F32); - } - psImageInit(out, 0.0); - - // - // Generate the spatial-order polynomial. The [i][j]-th element of - // the psImage polyValues will hold (x^i * y^j) for the stamp. - // - psImage *polyValues = GenSpatialOrder(spatialOrder, x, y); - - // XXX: switch (i, j) so they correspond to (x, y). - if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { - for (psS32 k = 0 ; k < nBF ; k++) { - psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; - psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; - psF64 polyVal = polyValues->data.F64[yOrder][xOrder]; - - // XXX: Verify that this is correct. - for (psS32 i = -kernelSize ; i <= kernelSize ; i++) { - for (psS32 j = -kernelSize ; j <= kernelSize ; j++) { - psImage *preCalc = (psImage *) kernels->preCalc->data[k]; - out->data.F32[i+kernelSize][j+kernelSize]+= - solution->data.F64[k] * - preCalc->data.F32[i+kernelSize][j+kernelSize] * - polyVal; - } - } - } - } else if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { - for (psS32 k = 0 ; k < nBF ; k++) { - // XXX: Why don't we have compilation warnings on type here (if - // we remove the (psS32) cast)? - psS32 u = (psS32) kernels->u->data.F32[k]; - psS32 v = (psS32) kernels->v->data.F32[k]; - psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; - psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; - // XXX: Verify that this is correct. - - out->data.F32[kernelSize - v][kernelSize - u]+= - solution->data.F64[k] * polyValues->data.F64[yOrder][xOrder]; - } - } - psFree(polyValues); - - psTrace("psModules.imcombine", 3, - "Exiting pmSubtractionKernelImage()\n"); - return(out); -} Index: unk/psModules/src/imcombine/pmImageSubtract.h =================================================================== --- /trunk/psModules/src/imcombine/pmImageSubtract.h (revision 13339) +++ (revision ) @@ -1,120 +1,0 @@ -/* @file ImageSubtract.h - * - * This file will contain code which creates a set of kernel basis - * functions, solves for their solution, and applies them to an image. - * - * @author Paul Price, IfA (original prototype) - * @author GLG, MHPCC - * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-30 21:12:56 $ - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PM_IMAGE_COMBINE_H -#define PM_IMAGE_COMBINE_H - -/// @addtogroup imcombine Image Combinations -/// @{ - -typedef enum { - PM_SUBTRACTION_KERNEL_POIS, ///< POIS kernel --- delta functions - PM_SUBTRACTION_KERNEL_ISIS ///< ISIS kernel --- gaussians modified by polynomials -} pmSubtractionKernelsType; - -typedef struct -{ - pmSubtractionKernelsType type; ///< Type ofKernels --- 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 *xOrder, *yOrder; ///< Spatial Polynomial order (for all) - int subIndex; ///< Index of kernel to be subtracted to maintain flux conservation - psArray *preCalc; ///< Array of images containing pre-calculated kernel (to - ///< accelerate ISIS; don't use for POIS) - psS32 size; ///< The halfsize of the kernel - psS32 spatialOrder; ///< The spatial order of the kernels -} -psSubtractionKernels; - -psSubtractionKernels *pmSubtractionKernelsAllocPOIS( - int size, - int SpatialOrder -); - -psSubtractionKernels *pmSubtractionKernelsAllocISIS( - const psVector *sigmas, - const psVector *orders, - int size, - int SpatialOrder -); - -typedef enum { - PM_STAMP_INIT, ///< Initial state - PM_STAMP_USED, ///< Use this stamp - PM_STAMP_REJECTED, ///< This stamp has been rejected - PM_STAMP_RECALC, ///< Having been reset, this stamp is to be recalculated - PM_STAMP_NONE ///< No stamp in this region -} pmStampStatus; - -typedef struct -{ - int x, y; ///< Position - int p_xSize; - int p_ySize; - int p_xMin; - int p_xMax; - int p_yMin; - int p_yMax; - psImage *matrix; ///< Associated matrix - psVector *vector; ///< Assoicated vector - pmStampStatus status; ///< Status ofstamp -} -pmStamp; - -psArray *pmSubtractionFindStamps( - psArray *stamps, ///< Output stamps, or NULL - const psImage *image, ///< Image for which to find stamps - const psImage *mask, ///< Mask - psU32 maskVal, ///< Value for mask - psF32 threshold, ///< Threshold for stamps in the image - psS32 xNum, ///< Number of stamps in x - psS32 yNum, ///< Number of stamps in y - psS32 border ///< Border around image to ignore (should be size of kernel) -); - -bool pmSubtractionCalculateEquation( - psArray *stamps, ///< The stamps for which to calculate the equation, - const psImage *reference, ///< Reference image - const psImage *input, ///< Input image - const psSubtractionKernels *kernels,///< The kernel basis functions - psS32 footprint ///< Half-size of region over which to calculate equation -); - - -psVector *pmSubtractionSolveEquation( - psVector *solution, ///< Solution vector, or NULL - const psArray *stamps ///< Array of stamps -); - -bool pmSubtractionRejectStamps( - psArray *stamps, ///< Array of stamps to check for rejection - psImage *mask, ///< Mask image - psU32 badStampMaskVal, ///< Value to use in mask for bad stamp - psS32 footprint, ///< Region to mask if stamp is bad - psF32 sigmaRej, ///< Number of RMS deviations above zero at which to reject - const psImage *refImage, ///< Reference image - const psImage *inImage, ///< Input image - const psVector *solution, ///< Solution vector - const psSubtractionKernels *kernels ///< Array of kernel parameters -); - -psImage *pmSubtractionKernelImage( - psImage *out, - const psVector *solution, - const psSubtractionKernels *kernels, - psF32 x, - psF32 y -); - -/// @} -#endif Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 13340) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 13340) @@ -0,0 +1,629 @@ +/** @file pmSubtraction.c + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-11 02:21:01 $ + * + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmHDU.h" // Required for pmFPA.h +#include "pmFPA.h" +#include "pmSubtractionStamps.h" + +#include "pmSubtraction.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Given a stamp coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j +static psImage *spatialPolyValues(int spatialOrder, // Maximum spatial polynomial order + float x, float y // Normalised position of interest, [-1,1] + ) +{ + assert(spatialOrder >= 0); + assert(x >= -1 && x <= 1); + assert(y >= -1 && y <= 1); + + psImage *polyValues = psImageAlloc(spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); // Matrix to return + psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, spatialOrder, spatialOrder); // Polynomial + + for (int j = 0; j <= spatialOrder; j++) { + for (int i = 0; i <= spatialOrder - j; i++) { + poly->coeff[i][j] = 1.0; + polyValues->data.F64[j][i] = polyValues->data.F64[j][i] = psPolynomial2DEval(poly, x, y); + poly->coeff[i][j] = 0.0; + } + } + psFree(poly); + + return polyValues; +} + +// Generate the convolved pixel value +static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions + int index, // Kernel basis function index + int x, int y, // Pixel around which to convolve + const psImage *image, // Image to convolve + const psImage *polyValues, // Spatial polynomial values + float (*weightFunc)(float value) // Function for weighting + ) +{ + switch (kernels->type) { + case PM_SUBTRACTION_KERNEL_POIS: { + // Convolution with a delta function is just the value specified by the offset + int u = kernels->u->data.S32[index]; // Offset in x + int v = kernels->v->data.S32[index]; // Offset in y + int xOrder = kernels->xOrder->data.S32[index]; // Polynomial order in x + int yOrder = kernels->yOrder->data.S32[index]; // Polynomial order in y + return weightFunc(polyValues->data.F64[yOrder][xOrder]) * image->data.F32[y + v][x + u]; + } + case PM_SUBTRACTION_KERNEL_ISIS: { + int xOrder = kernels->xOrder->data.S32[index]; // Polynomial order in x + int yOrder = kernels->yOrder->data.S32[index]; // Polynomial order in y + psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel + int size = kernels->size; // Kernel half-size + double sum = 0.0; // Accumulated sum from convolution + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += weightFunc(kernel->kernel[v][u]) * image->data.F64[y + v][x + u]; + } + } + return weightFunc(polyValues->data.F64[yOrder][xOrder]) * sum; + } + default: + psAbort("Should never get here."); + } + return NAN; +} + +// Weighting function for use with convolvePixel: no weighting applied, suitable for combining image pixels +static inline float imageWeighting(float value) +{ + return value; +} + +// Weighting function for use with convolvePixel: weighting suitable for combining variances +static inline float varianceWeighting(float value) +{ + return PS_SQR(value); +} + + +// Mark a pixel as blank in the image, mask and weight +static inline void markBlank(psImage *image, // Image to mark as blank + psImage *mask, // Mask to mark as blank (or NULL) + psImage *weight, // Weight map to mark as blank (or NULL) + int x, int y, // Coordinates to mark blank + psMaskType blank // Blank mask value + ) +{ + image->data.F32[y][x] = NAN; + if (mask) { + mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + } + if (weight) { + weight->data.F32[y][x] = NAN; + } + return; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input, + const psImage *weight, const pmSubtractionKernels *kernels, int footprint) +{ + PS_ASSERT_ARRAY_NON_NULL(stamps, false); + PS_ASSERT_IMAGE_NON_NULL(reference, false); + PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(input, false); + PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(reference, input, false); + PS_ASSERT_PTR_NON_NULL(kernels, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); + if (weight) { + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, input, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); + } + if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); + } + PS_ASSERT_INT_NONNEGATIVE(footprint, false); + + int kernelSize = kernels->size; // Half-size of kernel + PS_ASSERT_INT_LARGER_THAN(footprint, kernelSize, false); + + // Image size + int numCols = reference->numCols; + int numRows = reference->numRows; + + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numKernels = kernels->u->n; // Number of kernel basis functions + int numParams = numKernels + 1; // Total number of parameters to solve for: coefficient of each kernel + // basis function, and a constant background offset. + int bgIndex = numKernels; // Index in matrix for the background + psVector *convolutions = psVectorAlloc(numKernels, PS_TYPE_F64); // Convolutions for each kernel + + // We iterate over each stamp, allocate the matrix and vectors if + // necessary, and then calculate those matrix/vectors. + for (int i = 0; i < stamps->n; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + psTrace("psModules.imcombine", 6, "Inspecting stamp %d\n", i); + + if (stamp->status == PM_SUBTRACTION_STAMP_CALCULATE) { + psTrace("psModules.imcombine", 6, "Calculating for stamp %d\n", i); + psImage *stampMatrix = stamp->matrix; // Least-squares matrix for this stamp + psVector *stampVector = stamp->vector; // Least-squares vector for this stamp + + if (!stampMatrix) { + stamp->matrix = stampMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + } + assert(stampMatrix->type.type == PS_TYPE_F64); + assert(stampMatrix->numCols == numParams && stampMatrix->numRows == numParams); + psImageInit(stampMatrix, 0.0); + + if (!stampVector) { + stamp->vector = stampVector = psVectorAlloc(numParams, PS_TYPE_F64); + } + assert(stampVector->type.type == PS_TYPE_F64); + assert(stampVector->n == numParams); + psVectorInit(stampVector, 0.0); + + // Pre-evaluated spatial polynomial values + psImage *polyValues = spatialPolyValues(spatialOrder, + 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols, + 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows); + + for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { + for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { + float invNoise2 = 1.0 / (weight ? weight->data.F32[y][x] : + input->data.F32[y][x]); // Inverse square noise + + // Generate the convolutions + for (int i = 0; i < numKernels; i++) { + convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues, + imageWeighting); + } + + if (spatialOrder > 0) { + // Subtract a nominated convolution, to force flux conservation + int subIndex = kernels->subIndex; // Index of convolution to subtract + for (int i = 0; i < subIndex; i++) { + convolutions->data.F64[i] -= convolutions->data.F64[subIndex]; + } + for (int i = subIndex + 1; i < numKernels; i++) { + convolutions->data.F64[i] -= convolutions->data.F64[subIndex]; + } + } + + // Generate the least-squares matrix and vector + // Upper diagonal only + for (int i = 0; i < numKernels; i++) { + for (int j = i; j < numKernels; j++) { + stampMatrix->data.F64[i][j] += convolutions->data.F64[i] * + convolutions->data.F64[j] * invNoise2; + } + stampVector->data.F64[i] += input->data.F32[y][x] * convolutions->data.F64[i] * + invNoise2; + + // Background term + stampMatrix->data.F64[i][bgIndex] += convolutions->data.F64[i] * invNoise2; + } + // Background only terms + stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; + stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; + } + } + psFree(polyValues); + + // Fill in lower diagonal of symmetric matrix + for (int i = 0; i < numKernels; i++) { + for (int j = 0; j < i; j++) { + stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i]; + } + stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex]; + } + + stamp->status = PM_SUBTRACTION_STAMP_USED; + } + } + + psFree(convolutions); + + return true; +} + + +psVector *pmSubtractionSolveEquation(psVector *solution, const psArray *stamps) +{ + PS_ASSERT_ARRAY_NON_NULL(stamps, NULL); + + // Check inputs, while summing the stamp matrices and vectors + long numParams = -1; // Number of parameters + for (int i = 0; i < stamps->n; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + PS_ASSERT_PTR_NON_NULL(stamp, NULL); + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + if (numParams == -1) { + numParams = stamp->vector->n; + } else { + PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, NULL); + } + PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, NULL); + PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, NULL); + PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, NULL); + } + if (numParams == -1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found."); + return NULL; + } + + if (solution) { + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_SIZE(solution, numParams, NULL); + } + + // Accumulate the least-squares matricies and vectors + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix for all stamps + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector for all stamps + psVectorInit(sumVector, 0.0); + psImageInit(sumMatrix, 0.0); + for (int i = 0; i < stamps->n; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + } + } + + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); + psFree(sumMatrix); + if (!luMatrix) { + psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + + solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + if (!solution) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + return NULL; + } + + return solution; +} + + +int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage, + psImage *mask, psMaskType badStampMaskVal, const psVector *solution, + int footprint, float sigmaRej, const pmSubtractionKernels *kernels) +{ + PS_ASSERT_ARRAY_NON_NULL(stamps, -1); + PS_ASSERT_IMAGE_NON_EMPTY(refImage, -1); + PS_ASSERT_IMAGE_TYPE(refImage, PS_TYPE_F32, -1); + PS_ASSERT_IMAGE_NON_EMPTY(inImage, -1); + PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1); + PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -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); + PS_ASSERT_PTR_NON_NULL(kernels, -1); + PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, -1); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, -1); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, -1); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, -1); + if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + PS_ASSERT_ARRAY_NON_EMPTY(kernels->preCalc, -1); + PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, -1); + } + + int numKernels = solution->n - 1; // Number of kernel basis functions + + // Measure + psVector *deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); // Mean deviation for stamps + double totalSquareDev = 0.0; // Total square deviation from zero + int numStamps = 0; // Number of used stamps + { + int numCols = refImage->numCols, numRows = refImage->numRows; // Image dimensions + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics + + // Convolved image of stamp + psImage *convolvedStamp = psImageAlloc(2 * footprint + 1, 2 * footprint + 1, PS_TYPE_F32); + + for (int i = 0; i < stamps->n; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // The stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + int xStamp = stamp->x, yStamp = stamp->y; // Coordinates of stamp + psImageInit(convolvedStamp, 0.0); + + // Precalulate polynomial values + psImage *polyValues = spatialPolyValues(kernels->spatialOrder, + 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols, + 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows); + + for (int y = yStamp - footprint; y <= yStamp + footprint; y++) { + for (int x = xStamp - footprint; y <= xStamp + footprint; x++) { + for (int k = 0; k < numKernels; k++) { + convolvedStamp->data.F32[y - yStamp][x - xStamp] += + convolvePixel(kernels, k, x, y, refImage, polyValues, imageWeighting); + } + } + } + psFree(polyValues); + + psRegion region = psRegionSet(xStamp - footprint, xStamp + footprint, + yStamp - footprint, yStamp + footprint); // Region of interest + psImage *inStamp = psImageSubset(inImage, region); // Image of stamp + assert(inStamp->numCols == convolvedStamp->numCols && + inStamp->numRows == convolvedStamp->numRows); + (void)psBinaryOp(convolvedStamp, inStamp, "-", convolvedStamp); + (void)psBinaryOp(convolvedStamp, convolvedStamp, "/", inStamp); + psFree(inStamp); + (void)psBinaryOp(convolvedStamp, convolvedStamp, "*", convolvedStamp); + if (!psImageStats(stats, convolvedStamp, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, + "Unable to calculate statistics on normalised residual image of stamp."); + psFree(deviations); + psFree(convolvedStamp); + psFree(stats); + return -1; + } + deviations->data.F32[i] = stats->sampleMean; + totalSquareDev += PS_SQR(deviations->data.F32[i]); + numStamps++; + } + + psFree(convolvedStamp); + psFree(stats); + } + + float limit = sigmaRej * sqrt(totalSquareDev / (numStamps - 1)); // Limit on maximum deviation + + int numRejected = 0; + for (int i = 0; i < stamps->n; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + 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 + numRejected++; + for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { + for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { + mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal; + } + } + + // Set stamp for replacement + stamp->x = 0; + stamp->y = 0; + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + } + } + + psFree(deviations); + + return numRejected; +} + +psImage *pmSubtractionKernelImage(const psVector *solution, const pmSubtractionKernels *kernels, + float x, float y + ) +{ + PS_ASSERT_VECTOR_NON_NULL(solution, NULL); + PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); + PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); + if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, NULL); + PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, NULL); + } + + int numKernels = kernels->u->n; // Number of kernel basis functions + int size = kernels->size; // Size of kernel + + // We could simply sum the kernel basis functions, but convolving a delta function gives a more 'real' + // picture of what's actually happening in the convolution, and the speed difference can't be too large, + // since the image is small. + + psImage *image = psImageAlloc(2 * size + 1, 2 * size + 1, PS_TYPE_F32); // Output image + psImageInit(image, 0.0); + psImage *reference = psImageAlloc(4 * size + 1, 4 * size + 1, PS_TYPE_F32); // Reference image + psImageInit(reference, 0.0); + reference->data.F32[2 * size][2 * size] = 1.0; + + // Precalulate polynomial values + psImage *polyValues = spatialPolyValues(kernels->spatialOrder, x, y); + + for (int j = -size; j <= size; j++) { + for (int i = -size; i <= size; i++) { + for (int k = 0; k < numKernels; k++) { + image->data.F32[j + size][i + size] += solution->data.F64[k] * + convolvePixel(kernels, k, 2 * size + i, 2 * size + j, + reference, polyValues, imageWeighting); + } + } + } + + psFree(polyValues); + + return image; +} + +bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, + 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 (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) { + PS_ASSERT_IMAGE_NON_NULL(inWeight, false); + PS_ASSERT_IMAGE_TYPE(inWeight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(inWeight, inImage, false); + } + PS_ASSERT_PTR_NON_NULL(outImage, false); + if (*outImage) { + PS_ASSERT_IMAGE_NON_NULL(*outImage, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(*outImage, inImage, false); + PS_ASSERT_IMAGE_TYPE(*outImage, PS_TYPE_F32, false); + } + if (outMask && *outMask) { + PS_ASSERT_IMAGE_NON_NULL(*outMask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false); + PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGE_NON_NULL(inMask, false); + } + if (outWeight && *outWeight) { + PS_ASSERT_IMAGE_NON_NULL(*outWeight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(*outWeight, inImage, false); + PS_ASSERT_IMAGE_TYPE(*outWeight, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(inWeight, false); + } + PS_ASSERT_VECTOR_NON_NULL(solution, false); + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false); + PS_ASSERT_PTR_NON_NULL(kernels, false); + PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); + if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false); + PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, false); + } + + int numKernels = kernels->u->n; // Number of kernel basis functions + int size = kernels->size; // Half-size of kernel + float background = solution->data.F64[solution->n-1]; // The difference in background + int numCols = inImage->numCols, numRows = inImage->numRows; // Image dimensions + + psImage *convImage = *outImage; // Convolved image + if (!convImage) { + *outImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); + convImage = *outImage; + } + psImage *convMask = NULL; // Convolved mask image + if (outMask && inMask) { + if (!*outMask) { + *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + } + convMask = *outMask; + psImageInit(convMask, 0); + } + psImage *convWeight = NULL; // Convolved weight (variance) image + if (outWeight && inWeight) { + if (!*outWeight) { + *outWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight = *outWeight; + psImageInit(convWeight, 0.0); + } + + psImage *polyValues = NULL; // Pre-calculated polynomial values + int fullSize = 2 * size + 1; // Full size of kernel + + for (int y = size; y < inImage->numRows - size; y++) { + for (int x = size; x < inImage->numCols - size; x++) { + convImage->data.F32[y][x] = background; + + // Only generate polynomial values every kernel footprint, since we have already assumed + // (with the stamps) that it does not vary rapidly on this scale. + if (x % fullSize == size) { + psFree(polyValues); + polyValues = spatialPolyValues(kernels->spatialOrder, + 2.0 * (float)(x - numCols/2.0) / (float)numCols, + 2.0 * (float)(y - numRows/2.0) / (float)numRows); + } + + // Check and propagate the kernel footprint, if required + 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]; + } + } + 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; + } + } + + // Convolve the image + for (int i = 0; i < numKernels; i++) { + convImage->data.F32[y][x] += solution->data.F64[i] * + convolvePixel(kernels, i, x, y, inImage, polyValues, imageWeighting); + } + + // Convolve the weight (variance) map + if (inWeight) { + for (int i = 0; i < numKernels; i++) { + convWeight->data.F32[y][x] += PS_SQR(solution->data.F64[i]) * + convolvePixel(kernels, i, x, y, inWeight, polyValues, varianceWeighting); + } + } + } + } + psFree(polyValues); + + // Mark the non-convolved part as blank + for (int y = size; y < inImage->numRows - size; y++) { + for (int x = 0; x < size; x++) { + markBlank(convImage, convMask, convWeight, x, y, blank); + } + for (int x = numCols - size; x < numCols; x++) { + markBlank(convImage, convMask, convWeight, x, y, blank); + } + } + for (int y = 0; y < size; y++) { + for (int x = 0; x < numCols; x++) { + markBlank(convImage, convMask, convWeight, x, y, blank); + } + } + for (int y = numRows - size; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + markBlank(convImage, convMask, convWeight, x, y, blank); + } + } + + return true; +} Index: /trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.h (revision 13340) +++ /trunk/psModules/src/imcombine/pmSubtraction.h (revision 13340) @@ -0,0 +1,68 @@ +/* @file pmSubtraction.h + * + * PSF-matched image subtraction, based on the Alard & Lupton (1998) and Alard (2000) methods. + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-11 02:21:01 $ + * Copyright 2004-207 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_SUBTRACTION_H +#define PM_SUBTRACTION_H + +#include +#include "pmSubtractionKernels.h" + +/// @addtogroup imcombine Image Combinations +/// @{ + +/// Calculate the least-squares equation to match the image quality +bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, + const psImage *reference, ///< Reference image + const psImage *input, ///< Input image + const psImage *weight, ///< Weight image, or NULL + const pmSubtractionKernels *kernels, ///< Kernel parameters + int footprint ///< Half-size of region over which to calculate equation + ); + +/// Solve the least-squares equation to match the image quality +psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL + const psArray *stamps ///< Array of stamps + ); + +/// Reject stamps +int pmSubtractionRejectStamps(psArray *stamps, ///< Array of stamps to check for rejection + const psImage *refImage, ///< Reference image + psImage *inImage, ///< Input image + 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 + float sigmaRej, ///< Number of RMS deviations above zero at which to reject + const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// Generate an image of the convolution kernel +psImage *pmSubtractionKernelImage(const psVector *solution, ///< Solution vector + const pmSubtractionKernels *kernels, ///< Kernel parameters + float x, float y ///< Normalised position [-1,1] for which to generate image + ); + +/// Convolve image in preparation for subtraction +bool pmSubtractionConvolve(psImage **outImage, ///< Output image + psImage **outWeight, ///< Output weight map (or NULL) + psImage **outMask, ///< Output mask (or NULL) + const psImage *inImage, ///< Input image + const psImage *inWeight, ///< Input weight map (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 + const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// @} +#endif Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 13340) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 13340) @@ -0,0 +1,164 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmSubtractionKernels.h" + + +// Free function for pmSubtractionKernels +static void subtractionKernelsFree(pmSubtractionKernels *kernels) +{ + psFree(kernels->u); + psFree(kernels->v); + psFree(kernels->sigma); + psFree(kernels->xOrder); + psFree(kernels->yOrder); + psFree(kernels->preCalc); +} + +// Raise a number to an integer power +static inline float power(float value, // Value + int exp // Exponent + ) +{ + float result = value; // Result to return + for (int i = 2; i < exp; i++) { + result *= value; + } + return result; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, + int size, int spatialOrder) +{ + pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return + psMemSetDeallocator(kernels, (psFreeFunc)subtractionKernelsFree); + + kernels->type = type; + kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->sigma = NULL; + kernels->subIndex = 0; + kernels->preCalc = NULL; + kernels->size = size; + kernels->spatialOrder = spatialOrder; + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + + int num = PS_SQR(2 * size + 1) * (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, + size, spatialOrder); // The kernels + + // Generate a set of kernels for each (u,v) + for (int v = - size, index = 0; v <= size; v++) { + for (int u = - size; u <= size; u++) { + // 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->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + } + } + } + } + + kernels->subIndex = PS_SQR(size + 1); + assert(kernels->u->data.S32[kernels->subIndex] == 0 && + kernels->v->data.S32[kernels->subIndex] == 0 && + kernels->xOrder->data.S32[kernels->subIndex] == 0 && + kernels->yOrder->data.S32[kernels->subIndex] == 0); + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, + const psVector *sigmas, const psVector *orders) +{ + PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); + PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(sigmas, orders, NULL); + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + + int numGaussians = sigmas->n; // Number of Gaussians + + int num = 0; // Number of basis functions + for (int i = 0; i < numGaussians; i++) { + num += sigmas->data.F32[i] * orders->data.F32[i]; + } + num *= (spatialOrder + 1) * (spatialOrder + 2) / 2; + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, + size, spatialOrder); // The kernels + kernels->sigma = psVectorAlloc(num, PS_TYPE_F32); + kernels->preCalc = psArrayAlloc(num); + + // Set the kernel parameters + for (int i = 0, index = 0; i < numGaussians; i++) { + // Iterate over (u,v) order + for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { + for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) { + + // Set the pre-calculated kernel + psKernel *preCalc = psKernelAlloc(-size, size, -size, size); + for (int v = -size; v <= size; v++) { + for (int u = -size; u < size; u++) { + preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) * + expf(0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i])); + } + } + + // 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->sigma->data.F32[index] = sigmas->data.F32[i]; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + kernels->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc); + } + } + psFree(preCalc); // Drop reference + } + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, + const psVector *sigmas, const psVector *orders) +{ + switch (type) { + case PM_SUBTRACTION_KERNEL_POIS: + return pmSubtractionKernelsPOIS(size, spatialOrder); + case PM_SUBTRACTION_KERNEL_ISIS: + return pmSubtractionKernelsISIS(size, spatialOrder, sigmas, orders); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); + return NULL; + } +} Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 13340) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 13340) @@ -0,0 +1,54 @@ +#ifndef PM_SUBTRACTION_KERNELS_H +#define PM_SUBTRACTION_KERNELS_H + +#include + +/// Type of subtraction kernel +typedef enum { + PM_SUBTRACTION_KERNEL_POIS, ///< POIS kernel --- delta functions + PM_SUBTRACTION_KERNEL_ISIS ///< ISIS kernel --- gaussians modified by polynomials +} pmSubtractionKernelsType; + +/// Kernels specification +typedef struct { + 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 *xOrder, *yOrder; ///< Spatial Polynomial order (for all) + int subIndex; ///< Index of kernel to be subtracted (to maintain flux conservation) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) + int size; ///< The half-size of the kernel + int spatialOrder; ///< The spatial order of the kernels +} pmSubtractionKernels; + +/// General allocator for pmSubtractionKernels +/// +/// Unlike the functions for the specific kernel type, this function does not set up the basis functions, but +/// merely allocates space for their storage. +pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, ///< Number of basis functions + pmSubtractionKernelsType type, ///< Kernel type + int size, ///< Half-size of kernel + int spatialOrder ///< Order of spatial variations + ); + +/// Generate POIS kernels +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims) + int spatialOrder ///< Order of spatial variations + ); + +/// Generate ISIS kernels +pmSubtractionKernels *pmSubtractionKernelsISIS(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *sigmas, ///< Gaussian widths + const psVector *orders ///< Polynomial order of gaussians + ); + +/// Generate a kernel of a specified type +pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type + int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *sigmas, ///< Gaussian widths + const psVector *orders ///< Polynomial order of gaussians + ); + +#endif Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 13340) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 13340) @@ -0,0 +1,140 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmSubtractionStamps.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Free function for pmSubtractionStamp +void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free + ) +{ + psFree(stamp->matrix); + psFree(stamp->vector); +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +pmSubtractionStamp *pmSubtractionStampAlloc(pmSubtractionStampStatus status) +{ + pmSubtractionStamp *stamp = psAlloc(sizeof(pmSubtractionStamp)); // Stamp to return + psMemSetDeallocator(stamp, (psFreeFunc)subtractionStampFree); + + stamp->x = 0; + stamp->y = 0; + stamp->matrix = NULL; + stamp->vector = NULL; + stamp->status = status; + + return stamp; +} + + +psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *mask, + psMaskType maskVal, psMaskType used, float threshold, + float spacing, int footprint) +{ + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, 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(footprint, NULL); + PS_ASSERT_INT_LARGER_THAN(image->numCols, (2 * footprint), NULL); + PS_ASSERT_INT_LARGER_THAN(image->numRows, (2 * footprint), NULL); + + int xNumStamps = image->numCols / spacing; // Number of stamps in x dimension + int yNumStamps = image->numRows / spacing; // Number of stamps in y dimension + + if (stamps) { + PS_ASSERT_INT_EQUAL(stamps->n, xNumStamps * yNumStamps, NULL); + // Just in case, check for NULL stamps + for (int i = 0; i < xNumStamps * yNumStamps; i++) { + if (!stamps->data[i]) { + stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_REJECTED); + } + } + } else { + stamps = psArrayAlloc(xNumStamps * yNumStamps); + for (int i = 0; i < xNumStamps * yNumStamps ; i++) { + stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_INIT); + } + } + // Footprint of image + int numRows = image->numRows; + int numCols = image->numCols; + + for (int j = 0, index = 0; j < yNumStamps; j++) { + for (int i = 0; i < xNumStamps; i++, index++) { + pmSubtractionStamp *stamp = stamps->data[index]; // Stamp of interest + + // Only find a new stamp if we need to + if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || + stamp->status == PM_SUBTRACTION_STAMP_INIT) { + + // Find maximum non-masked value in the image section, + // but don't include a footprint around the edge + float fluxBest = threshold; // Flux of best stamp pixel + int xBest = 0, yBest = 0; // Coordinates of best stamp + + // Bounds of region to search for stamp + int yMin = footprint + j * (numCols - 2.0 * footprint) / yNumStamps; + int yMax = footprint + (j + 1) * (numCols - 2.0 * footprint) / yNumStamps - 1; + int xMin = footprint + i * (numRows - 2.0 * footprint) / xNumStamps; + int xMax = footprint + (i + 1) * (numRows - 2.0 * footprint) / xNumStamps - 1; + assert(yMax < image->numRows && xMax < image->numCols && yMin >= 0 && xMin >= 0); + + for (int y = yMin; y <= yMax ; y++) { + for (int x = xMin; x <= xMax ; x++) { + if (image->data.F32[y][x] > fluxBest) { + bool ok = true; + if (mask) { + // Check kernel footprint for bad pixels + if (mask->data.U8[y][x] & maskVal) { + ok = false; + } else { + for (int v = -footprint; v <= footprint && ok; v++) { + for (int u = -footprint; u <= footprint && ok; u++) { + if (mask->data.U8[y][x] & maskVal) { + ok = false; + } + } + } + } + } + + if (ok) { + fluxBest = image->data.F32[y][x]; + xBest = x; + yBest = y; + } + } + } + } + + stamp->x = xBest; + stamp->y = yBest; + stamp->status = fluxBest > threshold ? PM_SUBTRACTION_STAMP_CALCULATE : + PM_SUBTRACTION_STAMP_NONE; + } + } + } + + return stamps; +} Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 13340) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 13340) @@ -0,0 +1,36 @@ +#ifndef PM_SUBTRACTION_STAMPS_H +#define PM_SUBTRACTION_STAMPS_H + +#include + + +/// Status of stamp +typedef enum { + PM_SUBTRACTION_STAMP_INIT, ///< Initial state + PM_SUBTRACTION_STAMP_USED, ///< Use this stamp + PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected + PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp + PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region +} pmSubtractionStampStatus; + +/// A stamp for image subtraction +typedef struct { + int x, y; ///< Position + psImage *matrix; ///< Associated matrix + psVector *vector; ///< Assoicated vector + pmSubtractionStampStatus status; ///< Status ofstamp +} pmSubtractionStamp; + +/// Find stamps on an image +psArray *pmSubtractionFindStamps(psArray *stamps, ///< Output stamps, or NULL + 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 + int footprint ///< Half-size of stamp footprint + ); + + +#endif Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 13339) +++ /trunk/psModules/src/psmodules.h (revision 13340) @@ -11,5 +11,5 @@ // XXX the following headers define constructs needed by the elements below -#include +#include #include #include @@ -67,5 +67,7 @@ // the following headers are from psModule:imcombine #include -//#include +#include +#include +#include #include