Index: /trunk/psModules/src/pmImageSubtract.c =================================================================== --- /trunk/psModules/src/pmImageSubtract.c (revision 4033) +++ /trunk/psModules/src/pmImageSubtract.c (revision 4033) @@ -0,0 +1,1165 @@ +/** @file ImageSubtract.c + * + * This file will ... + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-25 23:00:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include "pslib.h" +#include "psConstants.h" +#include "pmImageSubtract.h" + +// XXX: Put this is psConstants.h +#define PS_ASSERT_VECTOR_SIZE(VEC, SIZE, RVAL) \ +if (VEC->n != SIZE) { \ + psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ + "psVector %s has size %d, should be %d." \ + #VEC, VEC->n, SIZE); \ + return(RVAL); \ +} + +/******************************************************************************* +XXX: Store the size and spatialOrder in a private member of the +psSubtractionKernels struct. + ******************************************************************************/ +psSubtractionKernels *pmSubtractionKernelsAllocPOIS(int size, + int 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_ISIS; + 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->p_size = size; + tmp->p_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 i = 0; i <= order; i++) { + psS32 j = order - i; + tmp->u->data.F32[ptr] = 0; + tmp->v->data.F32[ptr] = 0; + tmp->xOrder->data.F32[ptr] = i; + tmp->yOrder->data.F32[ptr] = j; + ptr++; + } + } + + // + // Iterate over (u,v). Generate a set of kernels foreach (u, v). + // + for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { + for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { + // 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^i * y^j such that (i+j) <= spatialOrder. + // + for (psS32 order = 0; order <= spatialOrder; order++) { + for (psS32 i = 0; i <= order; i++) { + psS32 j = order - i; + tmp->u->data.F32[ptr] = 0; + tmp->v->data.F32[ptr] = 0; + tmp->xOrder->data.F32[ptr] = i; + tmp->yOrder->data.F32[ptr] = j; + ptr++; + } + } + } + } + } + + return(tmp); +} + +/******************************************************************************* +XXX: Get the types correct. + +XXX: For the various spatial orders, verify that the preCalc data does not +vary with the (x^i * y^j) since that comes from the image coordinates when the +kernel is convolved with it. + +XXX: Code review this. + ******************************************************************************/ +psSubtractionKernels *pmSubtractionKernelsAllocISIS(const psVector *sigmas, + const psVector *orders, + int size, + int spatialOrder) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); + PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + // + // Calculate the number of basis functions (nBF), and others. + // + psS32 xKernelHalfSize = size; + psS32 yKernelHalfSize = size; + psS32 numSigmas = sigmas->n; + psS32 nBF = 0; + for (psS32 i = 0 ; i < numSigmas ; i++) { + nBF+= ((orders->n + 1) * (orders->n + 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); + // XXX: There should be a macro or a psLib vectorDuplicate(). + tmp->sigma = psVectorAlloc(numSigmas, PS_TYPE_F32); + for (psS32 i = 0 ; i < numSigmas ; i++) { + tmp->sigma->data.F32[i] = sigmas->data.F32[i]; + } + tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); + tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32); + // XXX: Is this correct? + tmp->subIndex = 0; + + tmp->p_size = size; + tmp->p_spatialOrder = spatialOrder; + + tmp->preCalc = psArrayAlloc(nBF); + for (psS32 i = 0 ; i < nBF ; i++) { + tmp->preCalc->data[i] = NULL; + } + + // + // This loop creates + // + // + // + // XXX: must loop on the sigmas. + psS32 sigPtr = 0; + // XXX: Review ptr. + psS32 kernelPtr = 0; + psArray *tmpPreCalc = psArrayAlloc(spatialOrder * spatialOrder); + for (psS32 order = 0; order <= spatialOrder; order++) { + for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { + psS32 orderYTerm = order - orderXTerm; + // psS32 imgPtr = orderYTerm + (orderXTerm * spatialOrder); + psImage *currKernel = (psImage *) psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); + + for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { + for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { + // Scale the (u,v) coordinates in kernel space to [-1.0:1.0]. + psF32 uScaled = ((psF32) (u - xKernelHalfSize)) / ((psF32) (1 + (2 * xKernelHalfSize))); + psF32 vScaled = ((psF32) (v - yKernelHalfSize)) / ((psF32) (1 + (2 * yKernelHalfSize))); + // 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])); + currKernel->data.F32[u+xKernelHalfSize][v+yKernelHalfSize] = + exp(exponent) * pow(uScaled, orderXTerm) + pow(vScaled, orderYTerm); + } + } + tmpPreCalc->data[kernelPtr++] = (psPtr *) currKernel; + } + } + + + // + // sigPtr loops through each sigma value in the *sigmas argument. + psS32 ptr = 0; + for (psS32 sigPtr = 0 ; sigPtr < numSigmas ; 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 orderP = 0 ; orderP < orders->data.S32[sigPtr] ; orderP++) { + for (psS32 xOrderP = 0 ; xOrderP < orderP ; xOrderP++) { + psS32 yOrderP = orderP - xOrderP; + + // + // We loop through all spatial orders. The next two loops have no effect + // on the values of the basis kernels and therefore, maybe, we should not + // be storing the preCalc data for them. + // + // XXX: Take this loop outside. Create the preCalc in a separate array. + // Then simply create pointers from the tmp struct to those preCalcs. + // + for (psS32 order = 0; order <= spatialOrder; order++) { + for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { + // XXX: Should simply add a pointer to an existing image and increase + // the memory reference counter for equivalent spatial orders. + if (NULL == tmp->preCalc->data[ptr]) { + tmp->preCalc->data[ptr] = (psPtr *) psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); + } else { + // XXX: gen warning. + printf("ERROR: tmp->preCalc->data[ptr] != NULL\n"); + } + + 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; + psImage *currKernel = (psImage *) tmp->preCalc->data[ptr]; + + for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { + for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { + // Scale the (u,v) coordinates in kernel space to [-1.0:1.0]. + psF32 uScaled = ((psF32) (u - xKernelHalfSize)) / ((psF32) (1 + (2 * xKernelHalfSize))); + psF32 vScaled = ((psF32) (v - yKernelHalfSize)) / ((psF32) (1 + (2 * yKernelHalfSize))); + + // 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])); + currKernel->data.F32[u+xKernelHalfSize][v+yKernelHalfSize] = + exp(exponent) * pow(uScaled, orderXTerm) + pow(vScaled, orderYTerm); + } + } + + ptr++; + } + } + } + } + } + + return(tmp); +} + +// XXX: These should maybe be public, maybe private global, or private local. +void p_pmStampFree(pmStamp *stamp) +{ + psFree(stamp->matrix); + psFree(stamp->vector); + psFree(stamp); +} + +// XXX: These should maybe be public, maybe private global, or private local. +pmStamp *p_pmStampAlloc(pmStampStatus status) +{ + pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp)); + stamp->x = 0; + stamp->y = 0; + stamp->matrix = NULL; + stamp->vector = NULL; + stamp->status = status; + + psMemSetDeallocator(stamp, (psFreeFcn)p_pmStampFree); + + return(stamp); +} + +/******************************************************************************* +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. + ******************************************************************************/ +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) + ) +{ + 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); + + if (stamps != NULL) { + PS_ASSERT_INT_EQUAL(stamps->n, (xNum * yNum), NULL); + // XXX: 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 + // + psS32 num = 0; + for (psS32 i = 0; i < xNum; i++) { + for (psS32 j = 0; j < yNum; j++) { + 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 over pixels within a + // border of pixels from the image edge. + // + // XXX: must set these + psS32 numX = 0; + psS32 numY = 0; + for (psS32 x = border + i * (numRows - 2.0 * border) / numX; + x < border + (i + 1) * (numRows - 2.0 * border) / numX; x++) { + for (psS32 y = border + j * (numCols - 2.0 * border) / numY; + y < border + (j + 1) * (numCols - 2.0 * border) / numY; y++) { + + // Determine if this pixel is larger than the max, and unmasked. + if ((image->data.F32[x][y] > max) && + !(mask->data.U8[x][y] & mask->data.U8[x][y])) { + max = image->data.F32[x][y]; + 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[bestx][besty] > 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++; + } + } + 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. + ******************************************************************************/ +static psImage GenSpatialOrder(psS32 spatialOrder, + psF32 x, + psF32 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*= imageY; + } + xSum*= imageX; + } + + return(polyValues); +} + + + +/******************************************************************************* + +XXX: How is the spatial order factor calculated? Is it simply a x^iy^j power +eavaluated at the (x, y) center for the stamp? Why bother with the 2-D +polynomial in evaluating it? Finally, why is it scaled to the half-size of +the reference image? + ******************************************************************************/ +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 + ) +{ + 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); + } + PS_ASSERT_INT_NONNEGATIVE(footprint, false); + psS32 numHalfRows = reference->numRows; + psS32 numHalfCols = reference->numCols; + psS32 spatialOrder = kernels->p_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 + + // + // The [i][j]-th element of this structure will hold (x^i * y^j) for the stamp. + // + psImage *polyValues = psImageAlloc(spatialOrder+1, spatialOrder+1, PS_TYPE_F64); + + // + // 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]; + + if (stamp->status == PM_STAMP_RECALC) { + psImage *stampMatrix = stamp->matrix; + psVector *stampVector = stamp->vector; + + if (stampMatrix == NULL) { + stampMatrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); + } else { + PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_SIZE(stampMatrix, numSolveParams, numSolveParams, false); + } + PS_IMAGE_SET_F64(stampMatrix, 0.0); + + if (stampVector == NULL) { + stampVector = psVectorAlloc(numSolveParams, PS_TYPE_F64); + } else { + PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, false); + PS_ASSERT_VECTOR_SIZE(stampVector, numSolveParams, false); + } + PS_VECTOR_SET_F64(stampVector, 0.0); + + // + // 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->y - numHalfCols)) / ((psF64) numHalfCols), + ((psF64) (stamp->x - numHalfRows)) / ((psF64) numHalfRows)); + + // + // Iterate over all pixels surrounding this stamp. + // + // + for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) { + for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) { + psF32 invNoise2 = 1.0/reference->data.F32[x][y]; // The inverse of the noise, squared. + + // XXX: Think about a better decomposition of the computation + // based on whether the kernels are ISIS or POIS. + + if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + // + // 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. + // + 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 = kernels->u->data.F32[k2]; // Offset in x + psS32 v2 = kernels->v->data.F32[k2]; // Offset in y + psS32 i2 = kernels->xOrder->data.F32[k2]; // Polynomial order in x + psS32 j2 = kernels->yOrder->data.F32[k2]; // Polynomial order in y + + // + // 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) { + printf("XXX: put some warning message here (ISIS kernels not implemented).\n"); + return(false); + // XXX: HEY: code this + /* + for (psS32 k1 = 0; k1 < numKernels; k1++) { + psF32 conv1 = GeneralKernelConvolve(reference, kernels, k1); + + for (psS32 k2 = k1; k2 < numKernels; k2++) { + psF32 conv2 = GeneralKernelConvolve(reference, kernels, k2); + + stampMatrix->data.F64[k1][k2] += conv1 * conv2 * invNoise2; + + } + vector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2; + stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2; + } + */ + } else { + printf("XXX: put some warning message here (bad kernel->type).\n"); + return(false); + } // if-else on kernel type. + } // y-loop for all pixels. + } // x-loop for all pixels. + psFree(polyValues); + + // + // 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: understand this. + // XXX: Why aren't they using numSolveParams instead of numKernels? + // XXX: is this POIS specific? + // XXX: What is config->penalty? + // + #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; + } + } + return(true); +} + + + + +/******************************************************************************* + ******************************************************************************/ +psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL + const psArray *stamps ///< Array of stamps + ) +{ + PS_ASSERT_PTR_NON_NULL(stamps, NULL); + psS32 size = ((pmStamp *) stamps->data[0])->vector->n; + + if (solution != NULL) { + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_SIZE(solution, size, NULL); + } else { + solution = psVectorAlloc(size, PS_TYPE_F64); + } + + psImage *sumMatrix = psImageAlloc(size, size, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(size, PS_TYPE_F64); + PS_VECTOR_SET_F64(sumVector, 0.0); + PS_IMAGE_SET_F64(sumMatrix, 0.0); + + for (psS32 s = 0; s < stamps->n; s++) { + pmStamp *stamp = (pmStamp *) stamps->data[s]; + psImage *stampMatrix = stamp->matrix; + psVector *stampVector = stamp->vector; + + PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_SIZE(stampVector, size, NULL); + PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, NULL); + PS_ASSERT_IMAGE_SIZE(stampMatrix, size, size, NULL); + + if (stamp->status == PM_STAMP_USED) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stampMatrix); + (void)psBinaryOp(sumVector, sumVector, "+", stampVector); + } + } + + psVector *permutation = NULL; + psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); + solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation); + + psFree(sumMatrix); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + + return(solution); +} + + +/******************************************************************************* +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: Parameter checks? + +XXX: Static structure for polyValues? + ******************************************************************************/ +statis psF32 *ConvolvePixelPois(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels + psS32 col, + psS32 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->p_spatialOrder; + + if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { + psF32 conv = background; // Initial convolved value + psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position + psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position + + // XXX: Ensure you have the (x, y) order correct. + 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 + // XXX: x/y order correct? + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + psF32 polyVal = polyValues->data.F32[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); + + 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: Parameter checks? + +XXX: Static structure for polyValues? + ******************************************************************************/ +statis psF32 *ConvolvePixelIsis(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels + psS32 col, + psS32 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->p_spatialOrder; + + if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) { + psF32 conv = background; // Initial convolved value + psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position + psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position + + // XXX: Ensure you have the (x, y) order correct. + psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); + + // Iterate over the kernel basis functions + for (psS32 k = 0; k < nBF; k++) { + // XXX: What's the story with this? + #if 0 + // XXX: x/y order correct? + psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; + psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; + psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; + #else + + psF32 polyVal = 1.0; + #endif + + psImage *preCalc = (psImage *) kernels->preCalc->data[nBF]; + for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { + for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { + // XXX: Should I do something special for the k==0 kernel? + // The POIS code does. + + conv += solution->data.F64[k] * + input->data.F32[yy+y][xx+x] * + preCalc->data.F32[yy-kernelSize][xx-kernelSize] * + polyVal; + } + } + } + } + psFree(polyValues); + + return(conv); +} + + + + + + +/******************************************************************************* +ConvolveImage(input, mask, badStampMaskVal, solution, kernels): convolves an +arbitrary image with either an ISIS or POIS set of kernel basis functions. + +XXX: Function decomposition: create a function which convolves a single +function with a set of POIS kernel basis functions. Create another for ISIS +kernel basis functions. Then rewrite this. + ******************************************************************************/ +psImage *ConvolveImage(const psImage *input, + const psImage *mask, + psU32 badStampMaskVal, + const psVector *solution, + const psSubtractionKernels *kernels) +{ + 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); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); + psS32 nBF = kernels->u->n; + PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); + + psS32 numCols = input->numCols; + psS32 numRows = input->numRows; + psS32 kernelSize = kernels->p_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 { + printf("XXX: Generate WARNING: unknown kernel type\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; + } + } + + return convolved; +} + + +/******************************************************************************* + ******************************************************************************/ +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) +{ + 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. + PS_VECTOR_SET_F32(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); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); + psS32 nBF = kernels->u->n; + PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); + + psS32 kernelSize = kernels->p_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) { + // XXX: Verify proper x,y order + psRegion myReg = psRegionSet(y - ySize, y + ySize, x - xSize, x + xSize); + 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); + // XXX: Verify proper x,y order + myReg = psRegionSet(kernelSize, kernelSize + 2 * footprint, + kernelSize, kernelSize + 2 * footprint); + psImage *subStampTrim = psImageSubset((psImage *) subStamp, myReg); + psImage *maskStampTrim = psImageSubset((psImage *) maskStamp, myReg); + // Copy image to workaround bug 305 + // XXX: Modify. This bug is now fixed. + // psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32); + // psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8); + // (void)psImageStats(stats, tempImage, tempMask, badStampMaskVal); + (void)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); + + 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 + ) +{ + 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); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); + psS32 nBF = kernels->u->n; + PS_ASSERT_VECTOR_SIZE(solution, 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 = PS_SQRT_F32(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); + return(true); +} + +/******************************************************************************* + ******************************************************************************/ +psImage *pmSubtractionKernelImage(psImage *out, + const psVector *solution, + const psSubtractionKernels *kernels, + psF32 x, + psF32 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); + // XXX: Assert that all vectors and arrays in kernels are of the same length. + psS32 nBF = kernels->u->n; + psS32 spatialOrder = kernels->p_spatialOrder; + psS32 kernelSize = kernels->p_size; + psS32 xCenter; // The pixel location for the center of the kernel in out img. + psS32 yCenter; // The pixel location for the center of the kernel in out img. + psS32 numRows; + psS32 numCols; + psS32 imgSize; + + if (out != NULL) { + xCenter = (psS32) ( ((x+1.0)/2.0) * (psF32) (out->numCols)); + yCenter = (psS32) ( ((y+1.0)/2.0) * (psF32) (out->numRows)); + if ( ((xCenter - kernelSize) < 0) || + ((xCenter + kernelSize) > out->numCols) || + ((yCenter - kernelSize) < 0) || + ((yCenter + kernelSize) > out->numRows)) { + + printf("XXX: generate WARNING: out image is not large enough.\n"); + return(out); + } + } else { + // + // We calculate the minimize size image that can contain the kernel + // centered at the specified (x, y) position. + // + numRows = (psS32) ceilf( ((psF32) kernelSize) / ((x + 1.0)/2.0)); + numRows = PS_MAX(numRows, (((psF32) kernelSize) / (1.0 - ((x + 1.0)/2.0)))); + numCols = (psS32) ceilf( ((psF32) kernelSize) / ((x + 1.0)/2.0)); + numCols = PS_MAX(numRows, (((psF32) kernelSize) / (1.0 - ((x + 1.0)/2.0)))); + imgSize = PS_MAX(numRows, numCols); + + out = psImageAlloc(imgSize, imgSize, PS_TYPE_F32); + xCenter = (psS32) ( ((x+1.0)/2.0) * (psF32) (out->numCols)); + yCenter = (psS32) ( ((y+1.0)/2.0) * (psF32) (out->numRows)); + } + PS_IMAGE_SET_F32(out, 0.0); + + if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { + for (psS32 k = 0 ; k < nBF ; k++) { + // XXX: verify these loop bounds with the size of the preCalc images. + for (psS32 i = -kernelSize ; i <= kernelSize ; i++) { + for (psS32 j = -kernelSize ; j <= kernelSize ; j++) { + psImage *preCalc = (psImage *) kernels->preCalc->data[nBF]; + out->data.F32[yCenter + i][xCenter + j]+= solution->data.F64[nBF] * + preCalc->data.F32[i+kernelSize][j+kernelSize]; + } + } + } + } else if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { + // + // 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(psS32 spatialOrder, x, y); + + for (psS32 k = 0 ; k < nBF ; k++) { + // XXX: Why don't we have compilation warnings on type here? + psS32 u = kernels->u->data.F32[nBF]; + psS32 v = kernels->v->data.F32[nBF]; + psS32 xOrder = kernels->xOrder->data.F32[nBF]; + psS32 yOrder = kernels->yOrder->data.F32[nBF]; + // XXX: Verify that (u, v) correspond to (x, y) and that we have the row,col stuff correct. + // XXX: Verify the same for xOrder,yOrder. + // XXX: Verify that this is correct. + out->data.F32[yCenter - (psS32) v][xCenter - (psS32) u]+= solution->data.F64[nBF] * + polyValues->data.F64[yOrder][xOrder]; + } + psFree(polyValues); + } + + return(out); +} +// This code is