Index: trunk/psModules/src/imsubtract/Makefile.am =================================================================== --- trunk/psModules/src/imsubtract/Makefile.am (revision 6873) +++ trunk/psModules/src/imsubtract/Makefile.am (revision 7019) @@ -3,11 +3,14 @@ libpsmoduleimsubtract_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) libpsmoduleimsubtract_la_LDFLAGS = -release $(PACKAGE_VERSION) -libpsmoduleimsubtract_la_SOURCES = pmImageSubtract.c \ - pmSubtractBias.c -# pmSubtractSky.c +libpsmoduleimsubtract_la_SOURCES = + +# pmImageSubtract.c \ +# pmSubtractBias.c +# pmSubtractSky.c psmoduleincludedir = $(includedir) -psmoduleinclude_HEADERS = \ - pmImageSubtract.h \ - pmSubtractBias.h +psmoduleinclude_HEADERS = + +# pmImageSubtract.h \ +# pmSubtractBias.h # pmSubtractSky.h Index: trunk/psModules/src/imsubtract/pmImageSubtract.c =================================================================== --- trunk/psModules/src/imsubtract/pmImageSubtract.c (revision 6873) +++ (revision ) @@ -1,1407 +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.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-04 01:01:33 $ - * - * 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. - */ - -#include -#include -#include "pslib.h" -#include "psConstants.h" -#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("ImageSubtract.pmSubtractionKernelsAllocPOIS", 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->u->n = tmp->u->nalloc; - tmp->v->n = tmp->v->nalloc; - tmp->sigma = NULL; - tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->xOrder->n = tmp->xOrder->nalloc; - tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32); - tmp->yOrder->n = tmp->yOrder->nalloc; - 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("ImageSubtract.pmSubtractionKernelsAllocPOIS", 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("ImageSubtract.pmSubtractionKernelsAllocISIS", 3, - "Calling pmSubtractionKernelsAllocISIS(%d, %d, %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->u->n = tmp->u->nalloc; - tmp->v->n = tmp->v->nalloc; - 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); - tmp->u->n = tmp->u->nalloc; - tmp->v->n = tmp->v->nalloc; - tmp->sigma->n = tmp->sigma->nalloc; - tmp->xOrder->n = tmp->xOrder->nalloc; - tmp->yOrder->n = tmp->yOrder->nalloc; - tmp->preCalc->n = tmp->preCalc->nalloc; - - // - // 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); - PS_IMAGE_SET_F32(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("ImageSubtract.pmSubtractionKernelsAllocISIS", 3, - "Exiting pmSubtractionKernelsAllocISIS(%d, %d, %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("ImageSubtract.pmSubtractionFindStamps", 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); - stamps->n = stamps->nalloc; - 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("ImageSubtract.pmSubtractionFindStamps", 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("ImageSubtract.GenSpatialOrder", 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("ImageSubtract.GenSpatialOrder", 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("ImageSubtract.IsisKernelConvolve", 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("ImageSubtract.IsisKernelConvolve", 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("ImageSubtract.ConvolvePixelPois", 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("ImageSubtract.ConvolvePixelPois", 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("ImageSubtract.ConvolvePixelIsis", 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("ImageSubtract.ConvolvePixelIsis", 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("ImageSubtract.ConvolveImage", 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, 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("ImageSubtract.ConvolveImage", 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("ImageSubtract.pmSubtractionCalculateEquation", 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("pmSubtractionCalculateEquation", 6, "subCalcEqn(): stamp %d\n", s); - if (stamp->status == PM_STAMP_RECALC) { - psTrace("pmSubtractionCalculateEquation", 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); - } - PS_IMAGE_SET_F64(stampMatrix, 0.0); - - if (stampVector == NULL) { - stampVector = psVectorAlloc(numSolveParams, PS_TYPE_F64); - stampVector->n = stampVector->nalloc; - stamp->vector = stampVector; - } else { - PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, false); - PS_ASSERT_VECTOR_SIZE(stampVector, numSolveParams, false); - } - PS_VECTOR_SET_F64(stampVector, 0.0); - psTrace("pmSubtractionCalculateEquation", 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("pmSubtractionCalculateEquation", 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("pmSubtractionCalculateEquation", 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("pmSubtractionCalculateEquation", 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("ImageSubtract.pmSubtractionCalculateEquation", 3, - "Exiting pmSubtractionCalculateEquation()\n"); - return(true); -} - - - - -/******************************************************************************* - ******************************************************************************/ -psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL - const psArray *stamps ///< Array of stamps - ) -{ - psTrace("ImageSubtract.pmSubtractionSolveEquation", 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, size, NULL); - } else { - solution = psVectorAlloc(size, PS_TYPE_F64); - solution->n = solution->nalloc; - } - - // - // 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); - sumVector->n = sumVector->nalloc; - PS_VECTOR_SET_F64(sumVector, 0.0); - PS_IMAGE_SET_F64(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, 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("ImageSubtract.pmSubtractionSolveEquation", 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("ImageSubtract.CalculateDeviations", 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); - deviations->n = deviations->nalloc; - // 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); - 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, 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("ImageSubtract.CalculateDeviations", 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("ImageSubtract.pmSubtractionRejectStamps", 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, 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("ImageSubtract.pmSubtractionRejectStamps", 3, - "Exiting pmSubtractionRejectStamps()\n"); - return(true); -} - -/******************************************************************************* - ******************************************************************************/ -psImage *pmSubtractionKernelImage(psImage *out, - const psVector *solution, - const psSubtractionKernels *kernels, - psF32 x, - psF32 y - ) -{ - psTrace("ImageSubtract.pmSubtractionKernelImage", 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); - } - PS_IMAGE_SET_F32(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("ImageSubtract.pmSubtractionKernelImage", 3, - "Exiting pmSubtractionKernelImage()\n"); - return(out); -} Index: trunk/psModules/src/imsubtract/pmImageSubtract.h =================================================================== --- trunk/psModules/src/imsubtract/pmImageSubtract.h (revision 6873) +++ (revision ) @@ -1,127 +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.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-28 20:43:52 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#if !defined(PM_IMAGE_COMBINE_H) -#define PM_IMAGE_COMBINE_H - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include -#include "pslib.h" -#include "psConstants.h" - -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/imsubtract/pmSubtractBias.c =================================================================== --- trunk/psModules/src/imsubtract/pmSubtractBias.c (revision 6873) +++ (revision ) @@ -1,577 +1,0 @@ -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the -// one that was being worked on. -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** @file pmSubtractBias.c - * - * This file will contain a module which will subtract the detector bias - * in place from an input image. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-17 18:10:08 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include "pmSubtractBias.h" - -#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 -#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 - - -#define MAX(a,b) ((a) > (b) ? (a) : (b)) -#define MIN(a,b) ((a) < (b) ? (a) : (b)) - - -// XXX: put these in psConstants.h -void PS_POLY1D_PRINT(psPolynomial1D *poly) -{ - printf("-------------- PS_POLY1D_PRINT() --------------\n"); - printf("poly->nX is %d\n", poly->nX); - for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) { - printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]); - } -} - -void PS_PRINT_SPLINE(psSpline1D *mySpline) -{ - printf("-------------- PS_PRINT_SPLINE() --------------\n"); - printf("mySpline->n is %d\n", mySpline->n); - for (psS32 i = 0 ; i < mySpline->n ; i++) { - PS_POLY1D_PRINT(mySpline->spline[i]); - } - PS_VECTOR_PRINT_F32(mySpline->knots); -} - -#define PS_IMAGE_PRINT_F32_HIDEF(NAME) \ -printf("======== printing %s ========\n", #NAME); \ -for (int i = 0 ; i < (NAME)->numRows ; i++) { \ - for (int j = 0 ; j < (NAME)->numCols ; j++) { \ - printf("%.5f ", (NAME)->data.F32[i][j]); \ - } \ - printf("\n"); \ -}\ - - -void overscanOptionsFree(pmOverscanOptions *options) -{ - psFree(options->stat); - psFree(options->poly); - psFree(options->spline); -} - -pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat) -{ - pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); - psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree); - - // Inputs - opts->single = single; - opts->fitType = fitType; - opts->order = order; - opts->stat = psMemIncrRefCounter(stat); - - // Outputs - opts->poly = NULL; - opts->spline = NULL; - - return opts; -} - - -/****************************************************************************** -psSubtractFrame(): this routine will take as input a readout for the input -image and a readout for the bias image. The bias image is subtracted in -place from the input image. -*****************************************************************************/ -static bool SubtractFrame(pmReadout *in,// Input readout - const pmReadout *sub, // Readout to be subtracted from input - float scale // Scale to apply before subtracting - ) -{ - assert(in); - assert(sub); - - psImage *inImage = in->image; // The input image - psImage *inMask = in->mask; // The input mask - psImage *subImage = sub->image; // The image to be subtracted - psImage *subMask = sub->mask; // The mask for the subtraction image - - // Offsets of the cells - int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0"); - int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0"); - int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0"); - int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0"); - - if ((inImage->numCols + x0in - x0sub) > subImage->numCols) { - psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n"); - return false; - } - if ((inImage->numRows + y0in - y0sub) > subImage->numRows) { - psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n"); - return false; - } - - if (scale == 1.0) { - for (int i = 0; i < inImage->numRows; i++) { - for (int j = 0; j < inImage->numCols; j++) { - inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub]; - if (inMask && subMask) { - inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; - } - } - } - } else { - for (int i = 0; i < inImage->numRows; i++) { - for (int j = 0; j < inImage->numCols; j++) { - inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale; - if (inMask && subMask) { - inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; - } - } - } - } - - return true; -} - - -#if 0 -/****************************************************************************** -ImageSubtractScalar(): subtract a scalar from the input image. - -XXX: Use a psLib function for this. - -XXX: This should - *****************************************************************************/ -static psImage *ImageSubtractScalar(psImage *image, - psF32 scalar) -{ - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - image->data.F32[i][j]-= scalar; - } - } - return(image); -} -#endif - -/****************************************************************************** -GenNewStatOptions(): this routine will take as input the options member of the -stat data structure, determine if multiple options have been specified, issue -a warning message if so, and return the highest priority option (according to -the order of the if-statements in this code). The higher priority options are -listed lower in the code. - *****************************************************************************/ -static psStatsOptions GenNewStatOptions(const psStats *stat) -{ - psS32 numOptions = 0; - psStatsOptions opt = 0; - - if (stat->options & PS_STAT_ROBUST_MEDIAN) { - if (numOptions == 0) { - opt = PS_STAT_ROBUST_MEDIAN; - } - numOptions++; - } - - if (stat->options & PS_STAT_CLIPPED_MEAN) { - if (numOptions == 0) { - opt = PS_STAT_CLIPPED_MEAN; - } - numOptions++; - } - - if (stat->options & PS_STAT_SAMPLE_MEDIAN) { - if (numOptions == 0) { - opt = PS_STAT_SAMPLE_MEDIAN; - } - numOptions++; - } - - if (stat->options & PS_STAT_SAMPLE_MEAN) { - numOptions++; - opt = PS_STAT_SAMPLE_MEAN; - } - - - if (numOptions == 0) { - psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n"); - } - if (numOptions != 1) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.c: GenNewStatOptions(): Too many statistics options have been specified\n"); - } - return(opt); -} - - - -#if 0 -/****************************************************************************** -ScaleOverscanVector(): this routine takes as input an arbitrary vector, -creates a new vector of length n, and fills the new vector with the -interpolated values of the old vector. The type of interpolation is: - PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data. - PM_FIT_SPLINE: fit splines to the input vector data. -XXX: Doesn't it make more sense to do polynomial interpolation on a few -elements of the input vector, rather than fit a polynomial to the entire -vector? - *****************************************************************************/ -static psVector *ScaleOverscanVector(psVector *overscanVector, - psS32 n, - void *fitSpec, - pmFit fit) -{ - psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, - "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n); - // PS_VECTOR_PRINT_F32(overscanVector); - - if (NULL == overscanVector) { - return(overscanVector); - } - - // Allocate the new vector. - psVector *newVec = psVectorAlloc(n, PS_TYPE_F32); - - // - // If the new vector is the same size as the old, simply copy the data. - // - if (n == overscanVector->n) { - for (psS32 i = 0 ; i < n ; i++) { - newVec->data.F32[i] = overscanVector->data.F32[i]; - } - return(newVec); - } - psPolynomial1D *myPoly; - psSpline1D *mySpline; - psF32 x; - psS32 i; - if (fit == PM_FIT_POLYNOMIAL) { - // Fit a polynomial to the old overscan vector. - myPoly = (psPolynomial1D *) fitSpec; - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); - if (myPoly == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n"); - return(NULL); - } - - // For each element of the new vector, convert the x-ordinate to that - // of the old vector, use the fitted polynomial to determine the - // interpolated value at that point, and set the new vector. - for (i=0;in) / ((psF32) n); - newVec->data.F32[i] = psPolynomial1DEval(myPoly, x); - } - } else if (fit == PM_FIT_SPLINE) { - psS32 mustFreeSpline = 0; - // Fit a spline to the old overscan vector. - mySpline = (psSpline1D *) fitSpec; - // XXX: Does it make any sense to have a psSpline argument? - if (mySpline == NULL) { - mustFreeSpline = 1; - } - - // - // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL, - // splines endpoints will be from 0.0 to overscanVector->n-1. Must scale - // properly when doing the spline eval. - // - // mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); - mySpline = psVectorFitSpline1D(NULL, overscanVector); - if (mySpline == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n"); - return(NULL); - } - // PS_PRINT_SPLINE(mySpline); - - // For each element of the new vector, convert the x-ordinate to that - // of the old vector, use the fitted polynomial to determine the - // interpolated value at that point, and set the new vector. - for (i=0;in - 1] - x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n); - newVec->data.F32[i] = psSpline1DEval(mySpline, x); - } - if (mustFreeSpline ==1) { - psFree(mySpline); - } - // PS_VECTOR_PRINT_F32(newVec); - - - } else { - psError(PS_ERR_UNKNOWN, true, "unknown fit type. Returning NULL.\n"); - psFree(newVec); - return(NULL); - } - - psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, - "---- ScaleOverscanVector() exit ----\n"); - return(newVec); -} - -#endif - -// Produce an overscan vector from an array of pixels -static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options - const psArray *pixels, // Array of vectors containing the pixel values - psStats *myStats // Statistic to use in reducing the overscan - ) -{ - // Reduce the overscans - psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row - psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate - psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting - for (int i = 0; i < pixels->n; i++) { - psVector *values = pixels->data[i]; // Vector with overscan values - if (values->n > 0) { - mask->data.U8[i] = 0; - ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] - psVectorStats(myStats, values, NULL, NULL, 0); - double reducedVal = NAN; // Result of statistics - if (! p_psGetStatValue(myStats, &reducedVal)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result " - "of statistics on row %d.\n", i); - return NULL; - } - reduced->data.F32[i] = reducedVal; - } else if (overscanOpts->fitType == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the " - "image, and no fit is requested.\n"); - return NULL; - } else { - // We'll fit this one out - mask->data.U8[i] = 1; - } - } - - // Fit the overscan, if required - switch (overscanOpts->fitType) { - case PM_FIT_NONE: - // No fitting --- that's easy. - break; - case PM_FIT_POLY_ORD: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - psFree(reduced); - reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_POLY_CHEBY: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - psFree(reduced); - reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_SPLINE: - // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an - // input spline - overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); - psFree(reduced); - reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); - return NULL; - break; - } - - psFree(ordinate); - psFree(mask); - - return reduced; -} - - - -/****************************************************************************** -XXX: The SDRS does not specify type support. F32 is implemented here. - *****************************************************************************/ -pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, - const pmReadout *bias, const pmReadout *dark) -{ - psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, - "---- pmSubtractBias() begin ----\n"); - PS_ASSERT_READOUT_NON_NULL(in, NULL); - PS_ASSERT_READOUT_NON_EMPTY(in, NULL); - PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); - - psImage *image = in->image; // The input image - - // Overscan processing - if (overscanOpts) { - // Check for an unallowable pmFit. - if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD && - overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) { - psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", overscanOpts->fitType); - return(in); - } - - psList *overscans = in->bias; // List of the overscan images - - psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original - myStats->options = GenNewStatOptions(overscanOpts->stat); - - // Reduce all overscan pixels to a single value - if (overscanOpts->single) { - psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); - pixels->n = 0; - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - int index = pixels->n; // Index - pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); - // XXX Reimplement with memcpy - for (int i = 0; i < overscan->numRows; i++) { - for (int j = 0; j < overscan->numCols; j++) { - pixels->data.F32[index++] = overscan->data.F32[i][j]; - } - } - - } - psFree(iter); - - (void)psVectorStats(myStats, pixels, NULL, NULL, 0); - double reduced = NAN; // Result of statistics - if (! p_psGetStatValue(myStats, &reduced)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning input image.\n"); - return(in); - } - (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32)); - } else { - - // We do the regular overscan subtraction - - bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction - - if (readRows) { - // The read direction is rows - psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - psVector *values = psVectorAlloc(0, PS_TYPE_F32); - values->n = 0; - pixels->data[i] = values; - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - int diff = image->row0 - overscan->row0; // Offset between the two regions - for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) { - // i is row on overscan - // XXX Reimplement with memcpy - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numCols); - for (int j = 0; j < overscan->numCols; j++) { - values->data.F32[index++] = overscan->data.F32[i][j]; - } - values->n += overscan->numCols; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = overscanVector(overscanOpts, pixels, myStats); - psFree(pixels); - if (! reduced) { - return in; - } - - // Subtract row by row - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - image->data.F32[i][j] -= reduced->data.F32[i]; - } - } - psFree(reduced); - - } else { - // The read direction is columns - psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - psVector *values = psVectorAlloc(0, PS_TYPE_F32); - values->n = 0; - pixels->data[i] = values; - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - int diff = image->col0 - overscan->col0; // Offset between the two regions - for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) { - // i is column on overscan - // XXX Reimplement with memcpy - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numRows); - for (int j = 0; j < overscan->numRows; j++) { - values->data.F32[index++] = overscan->data.F32[i][j]; - } - values->n += overscan->numRows; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = overscanVector(overscanOpts, pixels, myStats); - psFree(pixels); - if (! reduced) { - return in; - } - - // Subtract column by column - for (int i = 0; i < image->numCols; i++) { - for (int j = 0; j < image->numRows; j++) { - image->data.F32[j][i] -= reduced->data.F32[i]; - } - } - psFree(reduced); - } - } - psFree(myStats); - } // End of overscan subtraction - - // Bias frame subtraction - if (bias) { - SubtractFrame(in, bias, 1.0); - } - - if (dark) { - // Get the scaling - float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); - float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); - SubtractFrame(in, dark, inTime/darkTime); - } - - return in; -} - - Index: trunk/psModules/src/imsubtract/pmSubtractBias.h =================================================================== --- trunk/psModules/src/imsubtract/pmSubtractBias.h (revision 6873) +++ (revision ) @@ -1,77 +1,0 @@ -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the -// one that was being worked on. -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** @file pmSubtractBias.h - * - * This file will contain a module which will subtract the detector bias - * in place from an input image. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-17 18:10:08 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#if !defined(PM_SUBTRACT_BIAS_H) -#define PM_SUBTRACT_BIAS_H - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include -#include "pslib.h" - -#include "pmFPA.h" - -typedef enum { - PM_OVERSCAN_NONE, ///< No overscan subtraction - PM_OVERSCAN_EDGE, ///< Subtract the statistic of pixels along the to-be-determined readout direction - PM_OVERSCAN_ROWS, ///< Subtract rows - PM_OVERSCAN_COLUMNS, ///< Subtract columns - PM_OVERSCAN_ALL ///< Subtract the statistic of all pixels in overscan region -} pmOverscanAxis; - -typedef enum { - PM_FIT_NONE, ///< No fit - PM_FIT_POLY_ORD, ///< Fit ordinary polynomial - PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial - PM_FIT_SPLINE ///< Fit cubic splines -} pmFit; - -typedef struct -{ - // Inputs - bool single; // Reduce all overscan regions to a single value? - pmFit fitType; // Type of fit to overscan - unsigned int order; // Order of polynomial, or number of spline pieces - psStats *stat; // Statistic to use when reducing the minor direction - // Outputs - psPolynomial1D *poly; // Result of polynomial fit - psSpline1D *spline; // Result of spline fit -} -pmOverscanOptions; - -pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat); - -pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, - const pmReadout *bias, const pmReadout *dark); - -#if 0 -pmReadout *pmSubtractBias(pmReadout *in, ///< The input pmReadout image - void *fitSpec, ///< A polynomial or spline, defining the fit type. - const psList *overscans, ///< A psList of overscan images - pmOverscanAxis overScanAxis, ///< Defines how overscans are applied - psStats *stat, ///< The statistic to be used in combining overscan data - int nBin, ///< The amount of binning to be done image pixels. - pmFit fit, ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE - const pmReadout *bias); ///< A possibly NULL bias pmReadout which is to be subtracted -#endif - -#endif Index: trunk/psModules/src/imsubtract/pmSubtractSky.c =================================================================== --- trunk/psModules/src/imsubtract/pmSubtractSky.c (revision 6873) +++ (revision ) @@ -1,741 +1,0 @@ -/** @file pmSubtractSky.c - * - * This file will contain a module which will create a model of the - * background sky and subtract that from the input image. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-04 01:01:33 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - * - * - */ - -#include -#include -#include "pslib.h" -#include "pmSubtractSky.h" - -// XXX: Get rid of the. Create pmUtils.h -psImage *p_psDetermineTrimmedImage( - pmReadout *in -); - -/****************************************************************************** -DetermineNumBits(data): This routine takes an enum psStatsOptions as an -argument and returns the number of non-zero bits. - -XXX: This code is duplicated in the ReadoutCombine file. - *****************************************************************************/ -static psS32 DetermineNumBits(psStatsOptions data) -{ - psTrace("SubtractSky.DetermineNumBits", 4, "Calling DetermineNumBits(0x%x)\n", data); - - psS32 i; - psU64 tmpData = data; - psS32 numBits = 0; - - for (i=0;i<(8 * sizeof(psStatsOptions));i++) { - if (0x0001 & tmpData) { - numBits++; - } - tmpData = tmpData >> 1; - } - - psTrace("SubtractSky.DetermineNumBits", 4, - "Calling DetermineNumBits(0x%x) -> %d\n", data, numBits); - return(numBits); -} - -/****************************************************************************** -getHighestPriorityStatOption(statOptions): this routine takes as input a -psStats->options with multiple options set and returns one with a single -option set according to the precedence set in the SDRS. - *****************************************************************************/ -static psU64 getHighestPriorityStatOption(psU64 statOptions) -{ - psTrace("SubtractSky.getHighestPriorityStatOption", 4, - "Calling getHighestPriorityStatOption(0x%x)\n", statOptions); - - if (statOptions & PS_STAT_SAMPLE_MEAN) { - return(PS_STAT_SAMPLE_MEAN); - } else if (statOptions & PS_STAT_SAMPLE_MEDIAN) { - return(PS_STAT_SAMPLE_MEDIAN); - } else if (statOptions & PS_STAT_CLIPPED_MEAN) { - return(PS_STAT_CLIPPED_MEAN); - } else if (statOptions & PS_STAT_FITTED_MEAN) { - return(PS_STAT_FITTED_MEAN); - } else if (statOptions & PS_STAT_ROBUST_MEDIAN) { - return(PS_STAT_ROBUST_MEDIAN); - } - psError(PS_ERR_UNKNOWN, true, "Unallowable option requested for statistically binning image pixels.\n"); - return(-1); - // XXX - //else if (statOptions & PS_STAT_ROBUST_MODE) { - // return(PS_STAT_ROBUST_MODE); - //} -} - -/****************************************************************************** -psImage *binImage(origImage, binFactor, statOptions): This routine takes an -input psImage and scales it smaller by a factor of binFactor. The statistic -used in combining input pixels is specified in statOptions. - -XXX: use static vectors for myStats, binVector and binMask. -XXX: I coded this before I was aware of a psLib reBin function. I don't -use this function in this module. I'm keeping it here in the event that -requirements change and we might need a custom reBin function. - *****************************************************************************/ -/* -static psImage *binImage(psImage *origImage, - int binFactor, - psStatsOptions statOptions) -{ - psTrace("SubtractSky.binImage", 4, "Calling binImage(%d)\n", binFactor); - - if (binFactor <= 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: binImage(): binFactor is %d\n", binFactor); - return(origImage); - } - if (binFactor == 1) { - return(origImage); - } - - psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32); - psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8); - psStats *myStats = psStatsAlloc(statOptions); - - for (psS32 row = 0; row < origImage->numRows ; row+=binFactor) { - for (psS32 col = 0; col < origImage->numCols ; col+=binFactor) { - psS32 count = 0; - for (psS32 binRow = 0; binRow <= binFactor ; binRow++) { - for (psS32 binCol = 0; binCol <= binFactor ; binCol++) { - if (((row + binRow) < origImage->numRows) && - ((col + binCol) < origImage->numCols)) { - binVector->data.F32[count] = - origImage->data.F32[row + binRow][col + binCol]; - binMask->data.U8[count] = 0; - } else { - binVector->data.F32[count] = 0.0; - binMask->data.U8[count] = 1; - } - count++; - } - } - psStats *rc1 = psVectorStats(myStats, binVector, NULL, binMask, 1); - if (rc1 == NULL) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - return(origImage); - } - psF64 statValue; - psBool rc = p_psGetStatValue(rc1, &statValue); - - if (rc == true) { - origImage->data.F32[row][col] = (psF32) statValue; - } else { - origImage->data.F32[row][col] = 0.0; - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(), binImage(): p_psGetStatValue() was FALSE\n"); - } - } - } - psFree(binVector); - psFree(binMask); - psFree(myStats); - - psTrace("SubtractSky.binImage", 4, "Exiting binImage(%d)\n", binFactor); - return(origImage); -} -*/ - -/****************************************************************************** -CalculatePolyTerms(xOrder, yOrder): this routine will calculate the number of -coefficients (or terms) in a 2-D polynomial of order (xOrder, yOrder). - -XXX: Use your brain and figure out the analytical expression. - -XXX: Why isn't it simply (xOrder+1) * (yOrder+1)? - *****************************************************************************/ -static psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder) -{ - psTrace("SubtractSky.CalculatePolyTerms", 4, - "Calling CalculatePolyTerms(%d, %d)\n", xOrder, yOrder); - - psS32 maxOrder = PS_MAX(xOrder, yOrder); - psS32 localPolyTerms = 0; - psS32 order = 0; - psS32 num=0; - - for (order=0;order<=maxOrder;order++) { - for (num=0;num<=order;num++) { - if (((order-num) <= xOrder) && (num <= yOrder)) { - localPolyTerms++; - } - } - } - psTrace("SubtractSky.CalculatePolyTerms", 4, - "Exiting CalculatePolyTerms(%d, %d) -> %d\n", xOrder, yOrder, localPolyTerms); - return(localPolyTerms); - - // return((xOrder+1) * (yOrder+1)); -} - -/****************************************************************************** -buildPolyTerms(): this routine computes a 2-D array polyTerms[][] that holds -terms for the polynomial that is used to model the sky background. We use -this array primarily for convenience in computations involving sky model -polynomials. It is defined as: - polyTerms[i][0] = the power to which X is raised in the i-th term of in an - poly-order sky background polynomial. - - polyTerms[i][1] = the power to which Y is raised in the i-th term of in an - poly-order sky background polynomial. - *****************************************************************************/ -static psS32 **buildPolyTerms(psS32 xOrder, psS32 yOrder) -{ - psTrace("SubtractSky.buildPolyTerms", 4, - "Calling buildPolyTerms(%d, %d)\n", xOrder, yOrder); - - psS32 i=0; - psS32 order = 0; - psS32 num=0; - psS32 localPolyTerms = CalculatePolyTerms(xOrder, yOrder); - psS32 maxOrder = PS_MAX(xOrder, yOrder); - - // Create the data structure which we hold the xy order of each coeff. - psS32 **polyTerms = (psS32 **) psAlloc(localPolyTerms * sizeof(psS32 *)); - for (i=0; i < localPolyTerms ; i++) { - polyTerms[i] = (psS32 *) psAlloc(2 * sizeof(psS32)); - } - - i=0; - // This code segment loops through each term i in the polynomial and - // calculates the power to which x/y are raised in that i-th term. - // We first do the 0-order terms, then the 1-order terms, etc. - for (order=0;order<=maxOrder;order++) { - for (num=0;num<=order;num++) { - if (((order-num) <= xOrder) && (num <= yOrder)) { - polyTerms[i][0] = order-num; - polyTerms[i][1] = num; - i++; - } - } - } - - if (psTraceGetLevel(".psModule.pmSubtractSky.buildPolyTerms") >= 10) { - for (i=0; i < localPolyTerms ; i++) { - printf("x^%d * y^%d\n", polyTerms[i][0], polyTerms[i][1]); - } - } - - psTrace("SubtractSky.buildPolyTerms", 4, - "Exiting buildPolyTerms(%d, %d)\n", xOrder, yOrder); - return(polyTerms); -} - -/****************************************************************************** -This procedure calculates various combinations of powers of x and y and stores -them in the data structure p_psPolySums[][]. After it completes: - - p_psPolySums[i][j] == x^i * y^j - -XXX: Use a psImage for the p_psPolySums data structure? -XXX: p_psPolySums: should this be a global? Did you get the storage classifier - and name correct? -XXX: Use variable size arrays for p_psPolySums[][]. -XXX: Must initialize p_psPolySums[][]? - *****************************************************************************/ -#define PS_MAX_POLYNOMIAL_ORDER 20 - -psF64 p_psPolySums[PS_MAX_POLYNOMIAL_ORDER+1][PS_MAX_POLYNOMIAL_ORDER+1]; -static void buildSums(psF64 x, - psF64 y, - psS32 xOrder, - psS32 yOrder) -{ - psTrace("SubtractSky.buildPolyTerms", 4, - "Calling buildPolyTerms(%d, %d)\n", xOrder, yOrder); - - psS32 i = 0; - psS32 j = 0; - psF64 xSum = 0.0; - psF64 ySum = 0.0; - - xSum = 1.0; - ySum = 1.0; - for(i=0;i<=xOrder;i++) { - ySum = xSum; - for(j=0;j<=yOrder;j++) { - p_psPolySums[i][j] = ySum; - ySum*= y; - } - xSum*= x; - } - psTrace("SubtractSky.buildPolyTerms", 4, - "Exiting buildPolyTerms(%d, %d)\n", xOrder, yOrder); -} - -/****************************************************************************** -ImageFitPolynomial(myPoly, dataImage, maskImage): this private routine takes -an input image along with a mask and fits a polynomial to it. The degree of -the polynomial is specified by input parameter myPoly, and need not be -symmetrical in orders of X and Y. The polynomial must be type -PS_POLYNOMIAL_ORD. If there are not enough rows or columns in the input image -for the order of the polynomial, then that order is reduced. The algorithm -used in this routine is based on that of the pilot project ADD, but is not -documented anywhere. - -XXX: Different trace message facilities in use here. - *****************************************************************************/ -static psPolynomial2D *ImageFitPolynomial( - psPolynomial2D *myPoly, - psImage *dataImage, - psImage *maskImage) -{ - psTrace("SubtractSky.ImageFitPolynomial", 4, - "Calling ImageFitPolynomial()\n"); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_POLY_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_IMAGE_NON_NULL(dataImage, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(dataImage, NULL); - PS_ASSERT_IMAGE_TYPE(dataImage, PS_TYPE_F32, NULL); - PS_ASSERT_IMAGE_NON_NULL(maskImage, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(maskImage, NULL); - PS_ASSERT_IMAGE_TYPE(maskImage, PS_TYPE_U8, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(dataImage, maskImage, NULL); - psS32 oldPolyX = -1; - psS32 oldPolyY = -1; - - // The matrix equations become singular if there are more powers of X - // in myPoly then there are rows of the image. I think. Similarly for - // powers of Y and columns. So. Here we reduce the complexity of the - // polynomial if there are not enough rows/columns in the input image. - - if ((myPoly->nX + 1) > dataImage->numRows) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in x-dimension.\n"); - oldPolyX = myPoly->nX; - myPoly->nX = dataImage->numRows - 1; - } - if ((myPoly->nY + 1) > dataImage->numCols) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: ImageFitPolynomial(): Reducing polynomial complexity in y-dimension.\n"); - oldPolyY = myPoly->nY; - myPoly->nY = dataImage->numCols - 1; - } - psS32 i; - psS32 j; - psS32 x; - psS32 y; - psS32 aRow; - psS32 aCol; - psS32 **polyTerms = buildPolyTerms(myPoly->nX, myPoly->nY); - // We determine how many coefficients will be in the polynomial that we - // are fitting to this image. - psS32 localPolyTerms = CalculatePolyTerms(myPoly->nX, myPoly->nY); - psImage *A = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64); - psImage *Aout = psImageAlloc(localPolyTerms, localPolyTerms, PS_TYPE_F64); - psVector *B = psVectorAlloc(localPolyTerms, PS_TYPE_F64); - B->n = B->nalloc; - psVector *outPerm = NULL; - - // - // Initialize A matrix and B vector. - // - PS_IMAGE_SET_F64(A, 0.0); - PS_VECTOR_SET_F64(B, 0.0); - - // - // We build the A matrix and B vector. - // - for (x=0;xnumRows;x++) { - for (y=0;ynumCols;y++) { - if (maskImage->data.U8[x][y] == 0) { - buildSums((psF64) x, (psF64) y, myPoly->nX, myPoly->nY); - - /************************************************************ - This code dervies from equation (7) of the pilot ADD. However, - it is not exactly the same in that the order of the polynomial - may be different in X And Y. - - Equation (7) from the pilot ADD describes 16 linear equations. - The i-th equation is simply the partial derivative of the - sky background polynomial (1) w.r.t. to the i-th term in - that polynomial. The i-th equation is stored in row i of - matrix A[][] (matrix A[][] has origin (1,1), not (0,0)). To - compute A[i][j] we simply multiply the j-th term of the Sky - Background Polynomial (SBP) by the i-th term of SBP. - ************************************************************/ - for (aRow=0;aRowdata.F64[aRow][aCol]+= - (p_psPolySums[ polyTerms[aCol][0] ][ polyTerms[aCol][1] ] * - p_psPolySums[ polyTerms[aRow][0] ][ polyTerms[aRow][1] ]); - } - } - // Build the B[] vector, which is the right-hand side of (7). - for (i=0;idata.F64[i]+= dataImage->data.F32[x][y] * - p_psPolySums[ polyTerms[i][0] ][ polyTerms[i][1] ]; - } - } - } - } - - if (psTraceGetLevel(".psModule.pmSubtractSky.ImageFitPolynomial") >= 8) { - for (aRow=0;aRowdata.F64[aRow][aCol]); - } - } - - for (i=0;i<=localPolyTerms;i++) { - printf("B[%d] is %f\n", i, B->data.F64[i]); - } - } - - // - // Solve the matrix equations for the polynomial coefficients C. - // XXX: How do we know if these matrix operations were successful? - // - Aout = psMatrixLUD(Aout, &outPerm, A); - PS_ASSERT_IMAGE_NON_NULL(Aout, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(Aout, NULL); - psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64); - C->n = C->nalloc; - psMatrixLUSolve(C, Aout, B, outPerm); - - // - // Set the appropriate coefficients in the myPoly structure. - // - for (i=0;icoeff[ polyTerms[i][0] ][ polyTerms[i][1] ] = C->data.F64[i]; - psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6, - "myPoly->coeff[%d][%d] is %f\n", polyTerms[i][0], polyTerms[i][1], myPoly->coeff[ polyTerms[i][0] ][ polyTerms[i][1] ]); - } - - // - // Free data structures that were allocated in this module. - // - for (i=0;inX = oldPolyX; - for (i=oldPolyX ; i < (1 + myPoly->nX) ; i++) { - for (j=0;j<(1 + myPoly->nY) ; j++) { - myPoly->coeff[i][j] = 0.0; - } - } - } - if (oldPolyY != -1) { - myPoly->nY = oldPolyY; - for (i=0 ; i < (1 + myPoly->nX) ; i++) { - for (j=oldPolyY;j < (1 + myPoly->nY) ; j++) { - myPoly->coeff[i][j] = 0.0; - } - } - } - - psTrace("SubtractSky.ImageFitPolynomial", 4, - "Exiting ImageFitPolynomial()\n"); - // psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 4, - // "---- ImageFitPolynomial() end successfully ----\n"); - return(myPoly); -} - - -/****************************************************************************** -pmReadout pmSubtractSky(): - -XXX: use static vectors for myStats, and the binned image - -XXX: The SDR is silent about types. PS_TYPE_F32 is implemented here. - -XXX: Sync the psTrace message facilities. - *****************************************************************************/ -pmReadout *pmSubtractSky(pmReadout *in, - void *fitSpec, - psFit fit, - psS32 binFactor, - psStats *stats, - psF32 clipSD) -{ - PS_ASSERT_READOUT_NON_NULL(in, NULL); - PS_ASSERT_READOUT_NON_EMPTY(in, NULL); - PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } - psTrace(".psModule.pmSubtractSky", 4, - "---- pmSubtractSky() begin ----\n"); - - if ((fit != PM_FIT_NONE) && - (fit != PM_FIT_POLYNOMIAL) && - (fit != PM_FIT_SPLINE)) { - psError(PS_ERR_UNKNOWN, true, "psFit is unallowable (%d). Returning in image.\n", fit); - return(in); - } - - psStatsOptions statOptions = 0; - - // - // Return the original input readout if the fit specs are poorly defined. - // No warning or error messages should be generated. - // - if ((fitSpec == NULL) || - ((fit == PM_FIT_NONE) || (fit == PM_FIT_SPLINE))) { - // psLogMsg(__func__, PS_LOG_WARN, "Fit specs are poorly defined. Returning in image.\n"); - return(in); - } - - // - // Determine trimmed image from metadata. - // - - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - psImage *binnedImage = NULL; - psPolynomial2D *myPoly = NULL; - psImage *binnedMaskImage = NULL; - psU32 oldStatOptions = 0; - - // - // Determine which statistic to use when binning pixels, if any. - // - if (stats != NULL) { - statOptions = stats->options; - if (1 < DetermineNumBits(statOptions)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: Multiple statistical options have been requested.\n"); - statOptions = getHighestPriorityStatOption(statOptions); - if (statOptions == -1) { - psError(PS_ERR_UNKNOWN, true, "Not allowable stats->option was specified. Returning in image.\n"); - return(in); - } - // Save old input "stats" parameter. - oldStatOptions = stats->options; - stats->options = statOptions; - } - if (0 == DetermineNumBits(statOptions)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): no stats->options was requested\n"); - } - } - - // - // Generate required warning messages. - // - if (binFactor <= 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor); - } - if (stats == NULL) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): input parameter stats is NULL\n"); - } - - // - // Bin the input image according to input parameters. - // Create a new binned image mask. - // - if ((binFactor <= 1) || (stats == NULL) || (0 == DetermineNumBits(statOptions))) { - // No binning is required here. Simply create a copy of the image - // and a mask. - binnedImage = psImageCopy(binnedImage, trimmedImg, PS_TYPE_F32); - if (binnedImage == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageCopy() returned NULL. Returning in image.\n"); - return(in); - } - - if (in->mask != NULL) { - binnedMaskImage = psImageCopy(binnedMaskImage, in->mask, PS_TYPE_U8); - if (binnedMaskImage == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageCopy() returned NULL. Returning in image.\n"); - psFree(binnedImage); - return(in); - } - } else { - binnedMaskImage = psImageAlloc(binnedImage->numCols, - binnedImage->numRows, - PS_TYPE_U8); - PS_IMAGE_SET_U8(binnedMaskImage, 0); - } - } else { - binnedImage = psImageRebin(NULL, trimmedImg, in->mask, 0, binFactor, stats); - if (binnedImage == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageRebin() returned NULL. Returning in image.\n"); - return(in); - } - binnedMaskImage = psImageAlloc(binnedImage->numCols, - binnedImage->numRows, - PS_TYPE_U8); - PS_IMAGE_SET_U8(binnedMaskImage, 0); - } - psTrace(".psModule.pmSubtractSky", 4, - "binnedImage size is (%d, %d)\n", binnedImage->numRows, binnedImage->numCols); - - // - // Clip pixels that are outside the acceptable range. - // - if (clipSD <= 0.0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD); - } else { - // Determine the mean and standard deviation of the binned image. - psF64 binnedMean; - psF64 binnedStdev; - psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - psStats *rc = psImageStats(myStats, binnedImage, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); - return(in); - } - if (false == p_psGetStatValue(rc, &binnedMean)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine requested statistical operation. Returning in image.\n"); - return(in); - } - psTrace(".psModule.pmSubtractSky", 6, - "binned Mean is %f\n", binnedMean); - - myStats->options = PS_STAT_SAMPLE_STDEV; - rc = psImageStats(myStats, binnedImage, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); - return(in); - } - if (false == p_psGetStatValue(myStats, &binnedStdev)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine requested statistical operation. Returning in image.\n"); - return(in); - } - psFree(myStats); - psTrace(".psModule.pmSubtractSky", 6, - "binned StDev is %f\n", binnedStdev); - - // Clip all pixels which are more than clipSD sigmas from the mean. - psTrace(".psModule.pmSubtractSky", 6, - "clipSD is %f\n", clipSD); - - for (psS32 row = 0; row < binnedImage->numRows ; row++) { - for (psS32 col = 0; col < binnedImage->numCols ; col++) { - if (fabs(binnedImage->data.F32[row][col] - binnedMean) > - (clipSD * binnedStdev)) { - binnedMaskImage->data.U8[row][col] = 1; - } - } - } - } - - // - // Fit the polynomial to the binned image - // - if (fit == PM_FIT_POLYNOMIAL) { - myPoly = (psPolynomial2D *) fitSpec; - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_POLY_TYPE(myPoly, PS_POLYNOMIAL_ORD, NULL); - - myPoly = ImageFitPolynomial(myPoly, binnedImage, binnedMaskImage); - - if (myPoly != NULL) { - // Set the pixels in the binned image to that of the polynomial. - binnedImage = psImageEvalPolynomial(binnedImage, myPoly); - if (binnedImage == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageEvalPolynomial() returned NULL. Returning in image.\n"); - psFree(binnedMaskImage); - if (!((binFactor <= 1) || (stats == NULL))) { - psFree(binnedImage); - } - if (oldStatOptions != 0) { - stats->options = statOptions; - } - return(in); - } - } else { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractSky(): could not model sky with a polynomial.\n"); - psFree(binnedMaskImage); - if (!((binFactor <= 1) || (stats == NULL))) { - psFree(binnedImage); - } - if (oldStatOptions != 0) { - stats->options = statOptions; - } - return(in); - } - } else { - // We shouldn't get here since we check this above. - psError(PS_ERR_UNKNOWN, true, "Unallowable fit type. Returning in image.\n"); - psFree(binnedMaskImage); - if (!((binFactor <= 1) || (stats == NULL))) { - psFree(binnedImage); - } - if (oldStatOptions != 0) { - stats->options = statOptions; - } - return(in); - } - - // - //Subtract the polynomially fitted image from the original image - // - if (binFactor <= 1) { - // The binned image is the same size as the original image. - for (psS32 row = 0; row < trimmedImg->numRows ; row++) { - for (psS32 col = 0; col < trimmedImg->numCols ; col++) { - trimmedImg->data.F32[row][col]-= binnedImage->data.F32[row][col]; - } - } - } else { - for (psS32 row = 0; row < trimmedImg->numRows ; row++) { - for (psS32 col = 0; col < trimmedImg->numCols ; col++) { - // We calculate the F32 value of the pixel coordinates in the - // binned image and then use a pixel interpolation routine to - // determine the value of the pixel at that location. - psF32 binRowF64 = ((psF32) row) / ((psF32) binFactor); - psF32 binColF64 = ((psF32) col) / ((psF32) binFactor); - - // We add 0.5 to the pixel locations since the pixel - // interpolation routine defines the location of pixel - // (i, j) as (i+0.5, j+0.5). - binRowF64+= 0.5; - binColF64+= 0.5; - - psF32 binPixel = (psF32) psImagePixelInterpolate( - binnedImage, binColF64, binRowF64, - NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - trimmedImg->data.F32[row][col]-= binPixel; - - psTrace(".psModule.pmSubtractSky", 8, - "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f\n", - row, col, binRowF64-0.5, binColF64-0.5, binPixel); - } - } - - } - psFree(binnedMaskImage); - psFree(binnedImage); - if (oldStatOptions != 0) { - stats->options = statOptions; - } - - psTrace(".psModule.pmSubtractSky", 4, - "---- pmSubtractSky() exit successfully ----\n"); - return(in); -} Index: trunk/psModules/src/imsubtract/pmSubtractSky.h =================================================================== --- trunk/psModules/src/imsubtract/pmSubtractSky.h (revision 6873) +++ (revision ) @@ -1,40 +1,0 @@ -/** @file pmSubtractSky.h - * - * This file will contain a module which will create a model of the - * background sky and subtract that from the input image. - * - * @author GLG, MHPCC - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-17 18:10:08 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#if !defined(PM_SUBTRACT_SKY_H) -#define PM_SUBTRACT_SKY_H - -#if HAVE_CONFIG_H -#include -#endif - -#include -#include -#include "pslib.h" -#include "pmFPA.h" - -// XXX: this is pmFit in pmSubtractBias.c, named psFit here. -typedef enum { - PM_FIT_NONE, ///< No fit - PM_FIT_POLYNOMIAL, ///< Fit polynomial - PM_FIT_SPLINE ///< Fit cubic splines -} psFit; - -pmReadout *pmSubtractSky(pmReadout *in, - void *fitSpec, - psFit fit, - int binFactor, - psStats *stats, - float clipSD); -#endif