Index: trunk/psModules/src/imcombine/Makefile.am =================================================================== --- trunk/psModules/src/imcombine/Makefile.am (revision 14669) +++ trunk/psModules/src/imcombine/Makefile.am (revision 14671) @@ -11,5 +11,6 @@ pmSubtractionKernels.c \ pmSubtractionStamps.c \ - pmSubtractionMatch.c + pmSubtractionMatch.c \ + pmSubtractionParams.c pkginclude_HEADERS = \ @@ -20,5 +21,6 @@ pmSubtractionKernels.h \ pmSubtractionStamps.h \ - pmSubtractionMatch.h + pmSubtractionMatch.h \ + pmSubtractionParams.h CLEANFILES = *~ Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14669) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14671) @@ -40,50 +40,35 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Public functions +// Semi-public functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, - int size, int spatialOrder) -{ - pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return - psMemSetDeallocator(kernels, (psFreeFunc)subtractionKernelsFree); - - kernels->type = type; - kernels->description = NULL; - kernels->num = numBasisFunctions; - kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); - kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); - kernels->widths = NULL; - kernels->uStop = NULL; - kernels->vStop = NULL; - kernels->subIndex = 0; - kernels->preCalc = NULL; - kernels->size = size; - kernels->inner = 0; - kernels->spatialOrder = spatialOrder; - kernels->bgOrder = 0; - - return kernels; -} - -pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder) -{ - PS_ASSERT_INT_POSITIVE(size, NULL); - PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); - - int num = PS_SQR(2 * size + 1); // Number of basis functions - - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, - size, spatialOrder); // The kernels - psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", - size, spatialOrder, num); +bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, int start, int size) +{ + PS_ASSERT_PTR_NON_NULL(kernels, false); + PS_ASSERT_VECTOR_NON_NULL(kernels->u, false); + PS_ASSERT_VECTOR_NON_NULL(kernels->v, false); + PS_ASSERT_VECTOR_NON_NULL(kernels->widths, false); + PS_ASSERT_VECTOR_TYPE(kernels->u, PS_TYPE_S32, false); + PS_ASSERT_VECTOR_TYPE(kernels->v, PS_TYPE_S32, false); + PS_ASSERT_VECTOR_TYPE(kernels->widths, PS_TYPE_F32, false); + PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false); + PS_ASSERT_INT_NONNEGATIVE(start, false); + PS_ASSERT_INT_NONNEGATIVE(size, false); + + int numNew = PS_SQR(2 * size + 1); // Number of new kernel parameters to add + + // Ensure the sizes match + kernels->widths = psVectorRealloc(kernels->widths, start + numNew); + kernels->u = psVectorRealloc(kernels->u, start + numNew); + kernels->v = psVectorRealloc(kernels->v, start + numNew); + kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew); // Generate a set of kernels for each (u,v) - for (int v = - size, index = 0; v <= size; v++) { + for (int v = - size, index = start; v <= size; v++) { for (int u = - size; u <= size; u++, index++) { + kernels->widths->data.F32[index] = NAN; kernels->u->data.S32[index] = u; kernels->v->data.S32[index] = v; + kernels->preCalc->data[index] = NULL; psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); @@ -91,13 +76,9 @@ } - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0); - - return kernels; -} - -pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, - const psVector *fwhms, const psVector *orders) + return true; +} + +pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, + const psVector *fwhms, const psVector *orders) { PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); @@ -127,297 +108,5 @@ psFree(params); - kernels->widths = psVectorAlloc(num, PS_TYPE_F32); - kernels->preCalc = psArrayAlloc(num); - psKernel *subtract = NULL; // Kernel to subtract to maintain flux scaling - // Set the kernel parameters - for (int i = 0, index = 0; i < numGaussians; i++) { - float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma - float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian - - // Iterate over (u,v) order - for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { - for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { - - // Set the pre-calculated kernel - psKernel *preCalc = psKernelAlloc(-size, size, -size, size); - double sum = 0.0; // Normalisation - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * - expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma)); - } - } - if (index == 0) { - subtract = preCalc; - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - preCalc->kernel[v][u] /= sum; - } - } - } else if (uOrder % 2 == 0 && vOrder % 2 == 0) { - // Normalise sum of kernel component to unity for even functions - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u]; - } - } - } - - kernels->widths->data.F32[index] = fwhms->data.F32[i]; - kernels->u->data.S32[index] = uOrder; - kernels->v->data.S32[index] = vOrder; - kernels->preCalc->data[index] = preCalc; - - psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, - fwhms->data.F32[i], uOrder, vOrder); - } - } - } - - kernels->subIndex = -1; - - if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { - for (int i = 0; i < num; i++) { - psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest - psString kernelName = NULL; - psStringAppend(&kernelName, "kernel%03d.fits", i); - psFits *kernelFile = psFitsOpen(kernelName, "w"); - psFree(kernelName); - psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL); - psFitsClose(kernelFile); - } - } - - return kernels; -} - -/// Generate SPAM kernels -pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - int inner, ///< Inner radius to preserve unbinned - int binning ///< Kernel binning factor - ) -{ - PS_ASSERT_INT_POSITIVE(size, NULL); - PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); - PS_ASSERT_INT_NONNEGATIVE(inner, NULL); - PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); - PS_ASSERT_INT_POSITIVE(binning, NULL); - - // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region - int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region - int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region - assert(numOuter * binning + numInner == size); - int numTotal = numOuter + numInner; // Total number of summed pixels - - psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - - int num = PS_SQR(2 * numTotal + 1); // Number of basis functions - - psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); - - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, - size, spatialOrder); // The kernels - psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements", - size, inner, binning, spatialOrder, num); - - kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); - kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); - - psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element - psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element - locations->data.S32[numTotal] = 0; - widths->data.S32[numTotal] = 0; - for (int i = 1; i <= numInner; i++) { - locations->data.S32[numTotal + i] = i; - widths->data.S32[numTotal + i] = 0; - locations->data.S32[numTotal - i] = - i; - widths->data.S32[numTotal - i] = 0; - } - for (int i = numInner + 1; i <= numTotal; i++) { - locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] + - widths->data.S32[numTotal + i - 1] + 1; - widths->data.S32[numTotal + i] = binning - 1; - locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - binning; - widths->data.S32[numTotal - i] = binning - 1; - } - - if (psTraceGetLevel("psModules.imcombine") >= 10) { - for (int i = 0; i < 2 * numTotal + 1; i++) { - psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i], - locations->data.S32[i] + widths->data.S32[i]); - } - } - - // Set the kernel parameters - for (int i = - numTotal, index = 0; i <= numTotal; i++) { - int u = locations->data.S32[numTotal + i]; // Location of pixel - int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel - - for (int j = - numTotal; j <= numTotal; j++, index++) { - int v = locations->data.S32[numTotal + j]; // Location of pixel - int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel - - kernels->u->data.S32[index] = u; - kernels->v->data.S32[index] = v; - kernels->uStop->data.S32[index] = uStop; - kernels->vStop->data.S32[index] = vStop; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, - u, uStop, v, vStop); - } - } - - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0); - - psFree(locations); - psFree(widths); - - return kernels; -} - - -/// Generate FRIES kernels -pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - int inner ///< Inner radius to preserve unbinned - ) -{ - PS_ASSERT_INT_POSITIVE(size, NULL); - PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); - PS_ASSERT_INT_NONNEGATIVE(inner, NULL); - PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); - - int fibNum = 0; // Number of Fibonacci values - int fibLast = 1, fibTotal = 2; // Fibonacci sequence - while (fibTotal < size - inner) { - int temp = fibTotal; - fibTotal += fibLast; - fibLast = temp; - fibNum++; - } - - int numInner = inner; // Number of pixels in the inner region - int numOuter = fibNum; // Number of summed pixels in the outer region - int numTotal = numOuter + numInner; // Total number of summed pixels - - psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - - int num = PS_SQR(2 * numTotal + 1); // Number of basis functions - - psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); - - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, - size, spatialOrder); // The kernels - psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements", - size, inner, spatialOrder, num); - - kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); - kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); - - psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); - psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); - start->data.S32[numTotal] = 0; - stop->data.S32[numTotal] = 0; - for (int i = 1; i <= numInner; i++) { - start->data.S32[numTotal + i] = i; - stop->data.S32[numTotal + i] = i; - start->data.S32[numTotal - i] = -i; - stop->data.S32[numTotal - i] = -i; - } - for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal; - i++, fib = (temp = fib) + fibLast, fibLast = temp) { - start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1; - stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size); - start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i]; - stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i]; - } - - if (psTraceGetLevel("psModules.imcombine") >= 10) { - for (int i = 0; i < 2 * numTotal + 1; i++) { - psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]); - } - } - - // Set the kernel parameters - for (int i = - numTotal, index = 0; i <= numTotal; i++) { - int u = start->data.S32[numTotal + i]; // Location of pixel - int uStop = stop->data.S32[numTotal + i]; // Width of pixel - for (int j = - numTotal; j <= numTotal; j++, index++) { - int v = start->data.S32[numTotal + j]; // Location of pixel - int vStop = stop->data.S32[numTotal + j]; // Width of pixel - - kernels->u->data.S32[index] = u; - kernels->v->data.S32[index] = v; - kernels->uStop->data.S32[index] = uStop; - kernels->vStop->data.S32[index] = vStop; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, - u, uStop, v, vStop); - } - } - - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0); - - psFree(start); - psFree(stop); - - return kernels; -} - -// Grid United with Normal Kernel -pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, - const psVector *orders, int inner) -{ - PS_ASSERT_INT_POSITIVE(size, NULL); - PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); - PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); - PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(orders, NULL); - PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL); - PS_ASSERT_INT_NONNEGATIVE(inner, NULL); - PS_ASSERT_INT_LESS_THAN(inner, size, NULL); - - int numGaussians = fwhms->n; // Number of Gaussians - int numGaussianVars = 0; // Number of Gaussian variant functions in the kernel - psString params = NULL; // List of params - for (int i = 0; i < numGaussians; i++) { - int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian - numGaussianVars += (gaussOrder + 1) * (gaussOrder + 2) / 2; - psStringAppend(¶ms, "(%.2f,%d)", fwhms->data.F32[i], orders->data.S32[i]); - } - - int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements - - int num = numGaussianVars + numInner; // Total number of basis functions - - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK, - size, spatialOrder); // The kernels - psStringAppend(&kernels->description, "GUNK(%d,%d,%s,%d)", size, inner, params, spatialOrder); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "GUNK kernel: %d,%s,%d --> %d elements", - inner, params, spatialOrder, num); - psFree(params); - - kernels->widths = psVectorAlloc(numGaussianVars, PS_TYPE_F32); - kernels->preCalc = psArrayAlloc(numGaussianVars); - kernels->inner = numGaussianVars; - psKernel *subtract = NULL; // Kernel to subtract to maintain flux scaling - - // Set the Gaussian kernel parameters for (int i = 0, index = 0; i < numGaussians; i++) { float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma @@ -432,19 +121,13 @@ for (int u = -size; u <= size; u++) { sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * - expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(fwhms->data.F32[i])); + expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma)); } } - if (index == 0) { - subtract = preCalc; + + // Normalise sum of kernel component to unity for even functions + if (uOrder % 2 == 0 && vOrder % 2 == 0) { for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - preCalc->kernel[v][u] /= sum; - } - } - } else if (uOrder % 2 == 0 && vOrder % 2 == 0) { - // Normalise sum of kernel component to unity for even functions - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u]; + preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum; } } @@ -454,4 +137,7 @@ kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; + if (kernels->preCalc->data[index]) { + psFree(kernels->preCalc->data[index]); + } kernels->preCalc->data[index] = preCalc; @@ -462,15 +148,312 @@ } - // Generate a grid set of kernels for each (u,v) - for (int v = - inner, index = kernels->inner; v <= inner; v++) { - for (int u = - inner; u <= inner; u++, index++) { + kernels->subIndex = -1; + + if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { + for (int i = 0; i < num; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + psString kernelName = NULL; + psStringAppend(&kernelName, "kernel%03d.fits", i); + psFits *kernelFile = psFitsOpen(kernelName, "w"); + psFree(kernelName); + psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL); + psFitsClose(kernelFile); + } + } + + return kernels; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, + int size, int spatialOrder) +{ + pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return + psMemSetDeallocator(kernels, (psFreeFunc)subtractionKernelsFree); + + kernels->type = type; + kernels->description = NULL; + kernels->num = numBasisFunctions; + kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); + kernels->preCalc = psArrayAlloc(numBasisFunctions); + kernels->uStop = NULL; + kernels->vStop = NULL; + kernels->subIndex = 0; + kernels->size = size; + kernels->inner = 0; + kernels->spatialOrder = spatialOrder; + kernels->bgOrder = 0; + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + + int num = PS_SQR(2 * size + 1); // Number of basis functions + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, + size, spatialOrder); // The kernels + psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", + size, spatialOrder, num); + + if (!p_pmSubtractionKernelsAddGrid(kernels, 0, size)) { + psAbort("Should never get here."); + } + + kernels->subIndex = num / 2; + assert(kernels->u->data.S32[kernels->subIndex] == 0 && + kernels->v->data.S32[kernels->subIndex] == 0); + + return kernels; +} + + +pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, + const psVector *fwhms, const psVector *orders) +{ + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, + fwhms, orders); // Kernels + if (!kernels) { + return NULL; + } + + // Subtract a reference kernel from all the others to maintain flux scaling + psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from all others + for (int i = 1; i < kernels->num; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel to subtract + psBinaryOp(kernel->image, kernel->image, "-", subtract->image); + } + + return kernels; +} + +/// Generate SPAM kernels +pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner, ///< Inner radius to preserve unbinned + int binning ///< Kernel binning factor + ) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + PS_ASSERT_INT_POSITIVE(binning, NULL); + + // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region + int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region + int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region + assert(numOuter * binning + numInner == size); + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1); // Number of basis functions + + psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, + size, spatialOrder); // The kernels + psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements", + size, inner, binning, spatialOrder, num); + + kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + + psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element + psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element + locations->data.S32[numTotal] = 0; + widths->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + locations->data.S32[numTotal + i] = i; + widths->data.S32[numTotal + i] = 0; + locations->data.S32[numTotal - i] = - i; + widths->data.S32[numTotal - i] = 0; + } + for (int i = numInner + 1; i <= numTotal; i++) { + locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] + + widths->data.S32[numTotal + i - 1] + 1; + widths->data.S32[numTotal + i] = binning - 1; + locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - binning; + widths->data.S32[numTotal - i] = binning - 1; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i], + locations->data.S32[i] + widths->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = locations->data.S32[numTotal + i]; // Location of pixel + int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel + + for (int j = - numTotal; j <= numTotal; j++, index++) { + int v = locations->data.S32[numTotal + j]; // Location of pixel + int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel + kernels->u->data.S32[index] = u; kernels->v->data.S32[index] = v; - - psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); - } - } - - kernels->subIndex = numInner/2 + numGaussianVars; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, + u, uStop, v, vStop); + } + } + + kernels->subIndex = num / 2; + assert(kernels->u->data.S32[kernels->subIndex] == 0 && + kernels->v->data.S32[kernels->subIndex] == 0 && + kernels->uStop->data.S32[kernels->subIndex] == 0 && + kernels->vStop->data.S32[kernels->subIndex] == 0); + + psFree(locations); + psFree(widths); + + return kernels; +} + + +/// Generate FRIES kernels +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner ///< Inner radius to preserve unbinned + ) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + + int fibNum = 0; // Number of Fibonacci values + int fibLast = 1, fibTotal = 2; // Fibonacci sequence + while (fibTotal < size - inner) { + int temp = fibTotal; + fibTotal += fibLast; + fibLast = temp; + fibNum++; + } + + int numInner = inner; // Number of pixels in the inner region + int numOuter = fibNum; // Number of summed pixels in the outer region + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1); // Number of basis functions + + psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, + size, spatialOrder); // The kernels + psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements", + size, inner, spatialOrder, num); + + kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + + psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + start->data.S32[numTotal] = 0; + stop->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + start->data.S32[numTotal + i] = i; + stop->data.S32[numTotal + i] = i; + start->data.S32[numTotal - i] = -i; + stop->data.S32[numTotal - i] = -i; + } + for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal; + i++, fib = (temp = fib) + fibLast, fibLast = temp) { + start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1; + stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size); + start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i]; + stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i]; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = start->data.S32[numTotal + i]; // Location of pixel + int uStop = stop->data.S32[numTotal + i]; // Width of pixel + for (int j = - numTotal; j <= numTotal; j++, index++) { + int v = start->data.S32[numTotal + j]; // Location of pixel + int vStop = stop->data.S32[numTotal + j]; // Width of pixel + + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, + u, uStop, v, vStop); + } + } + + kernels->subIndex = num / 2; + assert(kernels->u->data.S32[kernels->subIndex] == 0 && + kernels->v->data.S32[kernels->subIndex] == 0 && + kernels->uStop->data.S32[kernels->subIndex] == 0 && + kernels->vStop->data.S32[kernels->subIndex] == 0); + + psFree(start); + psFree(stop); + + return kernels; +} + +// Grid United with Normal Kernel +pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, + const psVector *orders, int inner) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); + PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LESS_THAN(inner, size, NULL); + + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, + fwhms, orders); // Kernels + psStringPrepend(&kernels->description, "GUNK="); + psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); + + // Subtract unity from the kernels to maintain photometric flux scaling + for (int i = 0; i < kernels->num; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + kernel->kernel[0][0] -= 1.0; + } + + int numGaussians = kernels->num; // Number of ISIS kernels + int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements + + if (!p_pmSubtractionKernelsAddGrid(kernels, numGaussians, inner)) { + psAbort("Should never get here."); + } + + kernels->subIndex = numInner/2 + numGaussians; assert(kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0); @@ -518,6 +501,4 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements", size, inner, ringsOrder, spatialOrder, num); - - kernels->preCalc = psArrayAlloc(num); // Set the Gaussian kernel parameters Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14669) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14671) @@ -31,4 +31,10 @@ } pmSubtractionKernels; +/// Generate a delta-function grid for subtraction kernels (like the POIS kernel) +bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, ///< The subtraction kernels to append to + int start, ///< Index at which to start appending + int size ///< Half-size of the grid + ); + /// General allocator for pmSubtractionKernels /// @@ -44,4 +50,11 @@ pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims) int spatialOrder ///< Order of spatial variations + ); + +/// Generate ISIS kernels without the flux scaling built in +pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders ///< Polynomial order of gaussians ); Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 14669) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 14671) @@ -9,4 +9,5 @@ #include "pmHDU.h" #include "pmFPA.h" +#include "pmSubtractionParams.h" #include "pmSubtractionKernels.h" #include "pmSubtractionStamps.h" @@ -44,9 +45,77 @@ +static bool getStamps(psArray **stamps, // Stamps to read + const psArray *stampsData, // Stamp data from a file + const pmReadout *reference, // Reference readout + const pmReadout *input, // Input readout, or NULL to generate fake stamps + const psImage *subMask, // Mask for subtraction, or NULL + psImage *weight, // Weight map + const psRegion *region, // Region of interest, or NULL + float threshold, // Threshold for stamp finding + float stampSpacing, // Spacing between stamps + float targetWidth,// Target width for fake stamps + int size, // Kernel half-size + int footprint // Convolution footprint for stamps + ) +{ + + psImage *inImage = NULL; // Input image + if (input) { + inImage = input->image; + } + + if (stampsData) { + psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes + if (input) { + // We have x, y because the target is provided by the input image + xStamp = stampsData->data[0]; + yStamp = stampsData->data[1]; + } else { + // We have x, y and flux in order to generate a target + xStamp = stampsData->data[0]; + yStamp = stampsData->data[1]; + fluxStamp = stampsData->data[2]; + } + + psFree(*stamps); + *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region); + } else { + psTrace("psModules.imcombine", 3, "Finding stamps...\n"); + *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region, + threshold, stampSpacing); + } + if (!stamps) { + psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); + return false; + } + + memCheck(" find stamps"); + + if (!input && !pmSubtractionGenerateStamps(*stamps, targetWidth, footprint, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps."); + return false; + } + + psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); + if (!pmSubtractionExtractStamps(*stamps, reference->image, inImage, weight, footprint, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); + return false; + } + + memCheck(" extract stamps"); + + return true; +} + + + + + bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input, int footprint, float regionSize, float stampSpacing, float threshold, const char *stampsName, float targetWidth, pmSubtractionKernelsType type, - int size, int order, const psVector *isisWidths, const psVector *isisOrders, - int inner, int ringsOrder, int binning, int iter, float rej, psMaskType maskBad, + int size, int spatialOrder, const psVector *isisWidths, const psVector *isisOrders, + int inner, int ringsOrder, int binning, bool optimum, psVector *optFWHMs, + int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad, psMaskType maskBlank) { @@ -88,5 +157,5 @@ // We'll check kernel type when we allocate the kernels PS_ASSERT_INT_POSITIVE(size, false); - PS_ASSERT_INT_NONNEGATIVE(order, false); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, false); if (isisWidths || isisOrders) { PS_ASSERT_VECTOR_NON_NULL(isisWidths, false); @@ -99,4 +168,10 @@ PS_ASSERT_INT_NONNEGATIVE(ringsOrder, false); PS_ASSERT_INT_POSITIVE(binning, false); + if (optimum) { + PS_ASSERT_VECTOR_NON_NULL(optFWHMs, false); + PS_ASSERT_INT_NONNEGATIVE(optOrder, false); + PS_ASSERT_FLOAT_LARGER_THAN(optThreshold, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false); + } PS_ASSERT_INT_POSITIVE(iter, false); PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); @@ -136,4 +211,24 @@ } + // Read stamps from file + psArray *stampsData = NULL; // Stamps data read from file + if (stampsName && strlen(stampsName) > 0) { + psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName); + psVector *xStamp = NULL, *yStamp = NULL; // Stamp positions and fluxes + if (input) { + // We have x, y because the target is provided by the input image + stampsData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions + } else { + // We have x, y and flux in order to generate a target + stampsData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions + } + + // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions + xStamp = stampsData->data[0]; + yStamp = stampsData->data[1]; + psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + } + int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions @@ -145,10 +240,36 @@ memCheck("mask"); - // Kernel basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, isisWidths, isisOrders, - inner, binning, ringsOrder); + + psArray *stamps = NULL; // Stamps for matching PSF + + pmSubtractionKernels *kernels = NULL; // Kernel basis functions + if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { + if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, NULL, + threshold, stampSpacing, targetWidth, size, footprint)) { + goto ERROR; + } + kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, + stamps, footprint, optThreshold); + // XXX This is not quite the best thing to do here--- we'll find the same set of stamps again later, + // but since we may not be interested in the entire image on the first pass through, we have to blow + // the whole lot away. The alternative is to mark stamps that are outside the region of interest, and + // ignore them when generating sums and rejecting. + psFree(stamps); + stamps = NULL; + + if (!kernels) { + psErrorClear(); + psWarning("Unable to derive optimum ISIS kernel --- switching to default."); + } + } + + if (kernels == NULL) { + // Not an ISIS/GUNK kernel, or the optimum kernel search failed + kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, + inner, binning, ringsOrder); + } + psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", PS_DATA_UNKNOWN, "Subtraction kernels", kernels); - psArray *stamps = NULL; // Stamps for matching PSF psVector *solution = NULL; // Solution to match PSF @@ -186,51 +307,8 @@ psTrace("psModules.imcombine", 2, "Iteration %d...\n", k); - if (stampsName && strlen(stampsName) > 0) { - psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName); - iter = 1; // There is no iterating because we use all the stamps we have - psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes - if (input) { - // We have x, y because the target is provided by the input image - psArray *stampData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions - xStamp = stampData->data[0]; - yStamp = stampData->data[1]; - } else { - // We have x, y and flux in order to generate a target - psArray *stampData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions - xStamp = stampData->data[0]; - yStamp = stampData->data[1]; - fluxStamp = stampData->data[2]; - } - - // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions - psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - - stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region); - } else { - psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - stamps = pmSubtractionFindStamps(stamps, reference->image, subMask, region, - threshold, stampSpacing); - } - if (!stamps) { - psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); + if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, region, + threshold, stampSpacing, targetWidth, size, footprint)) { goto ERROR; } - - memCheck(" find stamps"); - - if (!input && !pmSubtractionGenerateStamps(stamps, targetWidth, footprint, kernels)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps."); - goto ERROR; - } - - psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionExtractStamps(stamps, reference->image, inImage, weight, - footprint, kernels)) { - psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); - goto ERROR; - } - - memCheck(" extract stamps"); psTrace("psModules.imcombine", 3, "Calculating equation...\n"); @@ -376,4 +454,6 @@ psFree(subMask); subMask = NULL; + psFree(stampsData); + stampsData = NULL; if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, kernels, maskBlank)) { @@ -395,4 +475,5 @@ psFree(stamps); psFree(solution); + psFree(stampsData); return false; } Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 14669) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 14671) @@ -28,4 +28,8 @@ int ringsOrder, ///< RINGS polynomial order int binning, ///< SPAM kernel binning + bool optimum, ///< Search for optimum ISIS kernel? + psVector *optFWHMs, ///< FWHMs for optimum search + int optOrder, ///< Maximum order for optimum search + float optThreshold, ///< Threshold for optimum search (0..1) // Operational parameters int iter, ///< Rejection iterations Index: trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 14671) +++ trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 14671) @@ -0,0 +1,416 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmSubtractionStamps.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Convolve the reference stamp by the kernel +static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved + const psKernel *kernel, // Kernel by which to convolve + int footprint // Size of area to be convolved + ) +{ + psKernel *convolution = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Result + psKernel *reference = stamp->reference; // Reference stamp, to be convolved + + // Range of kernel + int uMin = kernel->xMin; + int uMax = kernel->xMax; + int vMin = kernel->yMin; + int vMax = kernel->yMax; + + for (int y = -footprint; y <= footprint; y++) { + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, conv++) { + *conv = 0.0; + + int xStart = x + uMin; // Start index for convolution + for (int v = vMin; v <= vMax; v++) { + psF32 *ref = &reference->kernel[y + v][xStart]; // Dereference reference image + psF32 *krnl = &kernel->kernel[v][uMin]; // Dereference kernel in x + for (int u = uMin; u <= uMax; u++, ref++, krnl++) { + *conv += *ref * *krnl; + } + } + } + } + + return convolution; +} + + +// Accumulate cross-term sums for a stamp +static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2 + double *sumII, // Sum of I(x)^2/sigma(x)^2 + double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 + const pmSubtractionStamp *stamp, // Stamp with weight + const psKernel *input, // Input image, I(x) + const psKernel *convolution // Convolution, conv(x) + ) +{ + // Range of convolution + int xMin = convolution->xMin; + int xMax = convolution->xMax; + int yMin = convolution->yMin; + int yMax = convolution->yMax; + + psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + + for (int y = yMin; y <= yMax; y++) { + psF32 *in = &input->kernel[y][xMin]; // Dereference input + psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight + psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution + for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) { + double temp = *in / *wt; // Temporary product + *sumI += temp; + *sumII += *in * temp; + *sumIC += *conv * temp; + } + } + return; +} + +// Accumulate convolution sums for a stamp +static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2 + double *sumCC, // Sum of conv(x)^2/sigma(x)^2 + const pmSubtractionStamp *stamp, // Stamp with input and weight + const psKernel *convolution // Convolution, conv(x) + ) +{ + // Range of convolution + int xMin = convolution->xMin; + int xMax = convolution->xMax; + int yMin = convolution->yMin; + int yMax = convolution->yMax; + + psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + + for (int y = yMin; y <= yMax; y++) { + psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight + psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution + for (int x = xMin; x <= xMax; x++, wt++, conv++) { + double convNoise = *conv / *wt; // Temporary product + *sumC += convNoise; + *sumCC += *conv * convNoise; + } + } + return; +} + +static double accumulateChi2(psKernel *input, // Input stamp + psKernel *convolution, // Convolution, R(x)*k(u) + pmSubtractionStamp *stamp, // Stamp with weight + double coeff, // Coefficient of convolution + double bg // Background term + ) +{ + // Range of convolution + int xMin = convolution->xMin; + int xMax = convolution->xMax; + int yMin = convolution->yMin; + int yMax = convolution->yMax; + + double chi2 = 0.0; + psKernel *weight = stamp->weight; + + for (int y = yMin; y <= yMax; y++) { + psF32 *in = &input->kernel[y][xMin]; // Dereference input + psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight + psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution + for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) { + chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt; + } + } + + return chi2; +} + +// Return the initial value of chi^2 +static double initialChi2(psKernel *input, // Input stamp + const pmSubtractionStamp *stamp, // Stamp with weight + int footprint // Size of convolution + ) +{ + psKernel *weight = stamp->weight; // Weight map + psKernel *reference = stamp->reference; // Reference stamp + + double chi2 = 0.0; // Chi^2 + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *ref = &reference->kernel[y][-footprint]; // Derference reference + for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { + float diff = *in - *ref; // Temporary value + chi2 += PS_SQR(diff) / *wt; + } + } + + return chi2; +} + +// Subtract a convolution from the input +static void subtractConvolution(psKernel *input, // Input stamp + const pmSubtractionStamp *stamp, // Stamp with weight + const psKernel *convolution, // Convolution, R(x)*k(u) + float coeff, // Coefficient of subtraction + float bg // Background term + ) +{ + // Range of convolution + int xMin = convolution->xMin; + int xMax = convolution->xMax; + int yMin = convolution->yMin; + int yMax = convolution->yMax; + + psKernel *weight = stamp->weight; // Weight map + + for (int y = yMin; y <= yMax; y++) { + psF32 *in = &input->kernel[y][xMin]; // Dereference input + psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution + psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight + for (int x = xMin; x <= xMax; x++, in++, conv++, wt++) { + *in -= *conv * coeff + bg; + } + } + + return; +} + + +pmSubtractionKernels *pmSubtractionKernelsOptimumISIS(pmSubtractionKernelsType type, int size, int inner, + int spatialOrder, psVector *fwhms, int maxOrder, + const psArray *stamps, int footprint, float tolerance) +{ + if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Invalid kernel type: %x\n", type); + return NULL; + } + PS_ASSERT_INT_NONNEGATIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); + PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); + PS_ASSERT_INT_NONNEGATIVE(maxOrder, NULL); + PS_ASSERT_ARRAY_NON_NULL(stamps, NULL); + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, NULL); + + // Generate the kernels to test + int numGaussians = fwhms->n; // Number of Gaussians + int numKernels = numGaussians * (maxOrder + 1) * (maxOrder + 2) / 2; // Number of kernel components + psString params = NULL; // Parameter, for description + for (int i = 0; i < numGaussians; i++) { + psStringAppend(¶ms, "%.2f,", fwhms->data.F32[i]); + } + params[strlen(params) - 1] = '\0'; + + psVector *orders = psVectorAlloc(numGaussians, PS_TYPE_S32); // Polynomial orders + psVectorInit(orders, maxOrder); + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, + fwhms, orders); // Kernels + psFree(orders); + psFree(kernels->description); + kernels->description = NULL; + psStringAppend(&kernels->description, "OptISIS(%d,(%s),%d)", size, params, spatialOrder); + psFree(params); + + // Need to save the stamp inputs --- we're changing the values! + int numStamps = stamps->n; // Number of stamps + psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs + for (int i = 0; i < numStamps; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + psKernel *input = stamp->input; // Input image of interest + psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image + inputs->data[i] = psKernelAllocFromImage(copy, size + footprint + 1, size + footprint + 1); + psFree(copy); // Drop reference + } + + // Generate the convolutions + psArray *convolutions = psArrayAlloc(numKernels); // Convolutions per kernel component and stamp + for (int i = 0; i < numKernels; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Precalculated kernel + + psArray *stampConv = psArrayAlloc(numStamps); // Convolutions for each stamp + convolutions->data[i] = stampConv; + + for (int j = 0; j < numStamps; j++) { + stampConv->data[j] = convolveStamp(stamps->data[j], kernel, footprint); + } + } + + // Calculate sums invariant with the input image + double sum1 = 0.0; // sum of 1/sigma(x,y)^2 + for (int i = 0; i < numStamps; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + psKernel *weight = stamp->weight; // Weight map for stamp + + for (int v = -footprint; v <= footprint; v++) { + psF32 *wt = weight->kernel[v]; // Dereference weight map + for (int u = -footprint; u <= footprint; u++, wt++) { + sum1 += 1.0 / *wt; + } + } + } + psVector *sumC = psVectorAlloc(numKernels, PS_TYPE_F64); // sum of R(x)*k(u)/sigma(x)^2 + psVector *sumCC = psVectorAlloc(numKernels, PS_TYPE_F64); // sum of [R(x)*k(u)]^2/sigma(x)^2 + psVectorInit(sumC, 0.0); + psVectorInit(sumCC, 0.0); + for (int i = 0; i < numKernels; i++) { + psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp + for (int j = 0; j < numStamps; j++) { + accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j], + stampsConv->data[j]); + } + } + + // Initial chi^2 + double lastChi2 = 0.0; // Chi^2 from last iteration + int numPixels = numStamps * PS_SQR(2 * footprint + 1); // Number of pixels contributing to chi^2 + for (int i = 0; i < numStamps; i++) { + lastChi2 += initialChi2(inputs->data[i], stamps->data[i], footprint); + } + lastChi2 /= numPixels; + + // Rank the kernel components + psVector *ranking = psVectorAlloc(numKernels, PS_TYPE_S32); // Ranking of the kernel components + psVectorInit(ranking, -1); + int cutIndex = -1; // Index at which to cut off kernels + for (int iter = 0; iter < numKernels; iter++) { + int bestIndex = -1; // Index of best kernel component + double bestChi2 = INFINITY; // Value of best chi^2 + double bestCoeff = 0; // Value of best coefficient + double bestBG = 0; // Value of best background + + for (int i = 0; i < numKernels; i++) { + if (ranking->data.S32[i] >= 0) { + continue; + } + + double sumI = 0.0; // sum of I(x)/sigma(x)^2 + double sumII = 0.0; // sum of I(x)^2/sigma(x)^2 + double sumIC = 0.0; // sum of I(x)C(x)/sigma(x)^2 + + psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp + + for (int j = 0; j < numStamps; j++) { + accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], stampsConv->data[j]); + } + + double invDet = 1.0 / (sum1 * sumCC->data.F64[i] - PS_SQR(sumC->data.F64[i])); // Determinant^-1 + double coeff = invDet * (sum1 * sumIC - sumC->data.F64[i] * sumI); // Coefficient for kernel + double bg = invDet * (sumCC->data.F64[i] * sumI - sumC->data.F64[i] * sumIC); // Background + + double chi2 = 0.0; // Chi^2 + for (int j = 0; j < numStamps; j++) { + chi2 += accumulateChi2(inputs->data[j], stampsConv->data[j], stamps->data[j], coeff, bg); + } + + if (chi2 < bestChi2) { + bestIndex = i; + bestCoeff = coeff; + bestChi2 = chi2; + bestBG = bg; + } + + psTrace("psModules.imcombine", 8, "%d: %lf %lf %lf %lf %lf %lf\n", i, sum1, sumI, sumII, sumIC, + sumC->data.F64[i], sumCC->data.F64[i]); + psTrace("psModules.imcombine", 6, "%d: %lf %lf %lf\n", i, coeff, bg, chi2); + } + bestChi2 /= numPixels; + + // And the winner is.... + ranking->data.S32[bestIndex] = iter; + // Remove its contribution, and don't include it in the future. + psArray *stampsConv = convolutions->data[bestIndex]; // Convolutions per stamp + for (int j = 0; j < numStamps; j++) { + subtractConvolution(inputs->data[j], stamps->data[j], stampsConv->data[j], + bestCoeff, bestBG); + } + + double diff = lastChi2 - bestChi2; // Difference in chi^2 between iterations + + psTrace("psModules.imcombine", 3, "The winner of round %d is %d (%f,%d,%d): %lf (%lf) --> %lf\n", + iter, bestIndex, kernels->widths->data.F32[bestIndex], kernels->u->data.S32[bestIndex], + kernels->v->data.S32[bestIndex], bestCoeff, diff, bestChi2); + + if (fabsf(diff) < tolerance) { + cutIndex = iter; + break; + } + + lastChi2 = bestChi2; + } + psFree(inputs); + psFree(convolutions); + psFree(sumC); + psFree(sumCC); + + if (cutIndex < 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to converge to tolerance %g\n", tolerance); + psFree(ranking); + psFree(kernels); + return NULL; + } + + int newSize = cutIndex + 1; // Size of new kernel basis set + psTrace("psModules.imcombine", 2, "Accepting %d kernels.\n", newSize); + psVector *uNew = psVectorAlloc(newSize, PS_TYPE_S32); + psVector *vNew = psVectorAlloc(newSize, PS_TYPE_S32); + psVector *widthsNew = psVectorAlloc(newSize, PS_TYPE_F32); + psArray *preCalcNew = psArrayAlloc(newSize); + + for (int i = 0; i < numKernels; i++) { + int rank = ranking->data.S32[i]; // This kernel component's ranking + if (rank >= 0 && rank < newSize) { + uNew->data.S32[rank] = kernels->u->data.S32[i]; + vNew->data.S32[rank] = kernels->v->data.S32[i]; + widthsNew->data.F32[rank] = kernels->widths->data.F32[i]; + preCalcNew->data[rank] = psMemIncrRefCounter(kernels->preCalc->data[i]); + } + } + psFree(kernels->u); + psFree(kernels->v); + psFree(kernels->widths); + psFree(kernels->preCalc); + kernels->u = uNew; + kernels->v = vNew; + kernels->widths = widthsNew; + kernels->preCalc = preCalcNew; + kernels->num = newSize; + + psFree(ranking); + + // Maintain photometric scaling + if (type == PM_SUBTRACTION_KERNEL_ISIS) { + psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest + for (int i = 1; i < newSize; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + psBinaryOp(kernel->image, kernel->image, "-", subtract->image); + } + } else if (type == PM_SUBTRACTION_KERNEL_GUNK) { + psStringPrepend(&kernels->description, "GUNK="); + psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); + + for (int i = 0; i < newSize; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + kernel->kernel[0][0] -= 1.0; + } + + if (!p_pmSubtractionKernelsAddGrid(kernels, numGaussians, inner)) { + psAbort("Should never get here."); + } + + kernels->subIndex = PS_SQR(2 * inner + 1) / 2 + numGaussians; + assert(kernels->u->data.S32[kernels->subIndex] == 0 && + kernels->v->data.S32[kernels->subIndex] == 0); + } + + return kernels; +} Index: trunk/psModules/src/imcombine/pmSubtractionParams.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.h (revision 14671) +++ trunk/psModules/src/imcombine/pmSubtractionParams.h (revision 14671) @@ -0,0 +1,19 @@ +#ifndef PM_SUBTRACTION_PARAMS_H +#define PM_SUBTRACTION_PARAMS_H + +#include +#include "pmSubtractionKernels.h" + +/// Generate a set of optimum kernels for ISIS (or GUNK) +pmSubtractionKernels *pmSubtractionKernelsOptimumISIS(pmSubtractionKernelsType type, ///< Kernel type + int size, ///< Half-size of kernel + int inner, ///< Inner radius for GUNK + int spatialOrder, ///< Spatial polynomial order + psVector *fwhms, ///< Gaussian FWHMs to try + int maxOrder, ///< Maximum polynomial order + const psArray *stamps, ///< Stamps + int footprint, ///< Convolution footprint for stamps + float tolerance ///< Maximum difference in chi^2 + ); + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14669) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14671) @@ -239,5 +239,5 @@ bool pmSubtractionExtractStamps(psArray *stamps, psImage *reference, psImage *input, psImage *weight, - int footprint, const pmSubtractionKernels *kernels) + int footprint, int kernelSize) { PS_ASSERT_ARRAY_NON_NULL(stamps, false); @@ -254,7 +254,9 @@ PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); } + PS_ASSERT_INT_NONNEGATIVE(footprint, false); + PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); int numCols = reference->numCols, numRows = reference->numRows; // Size of images - int size = kernels->size + footprint; // Size of postage stamps + int size = kernelSize + footprint; // Size of postage stamps if (!weight) { @@ -303,6 +305,5 @@ -bool pmSubtractionGenerateStamps(psArray *stamps, float fwhm, int footprint, - const pmSubtractionKernels *kernels) +bool pmSubtractionGenerateStamps(psArray *stamps, float fwhm, int footprint, int kernelSize) { PS_ASSERT_ARRAY_NON_NULL(stamps, false); @@ -310,5 +311,5 @@ PS_ASSERT_INT_NONNEGATIVE(footprint, false); - int size = kernels->size + footprint; // Size of postage stamps + int size = kernelSize + footprint; // Size of postage stamps int num = stamps->n; // Number of stamps float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14669) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14671) @@ -51,5 +51,5 @@ psImage *weight, ///< Weight (variance) map int footprint, ///< Stamp footprint size - const pmSubtractionKernels *kernels ///< Kernel basis functions (for size) + int kernelSize ///< Kernel half-size ); @@ -58,5 +58,5 @@ float fwhm, ///< FWHM for each stamp int footprint, ///< Stamp footprint size - const pmSubtractionKernels *kernels ///< Kernel basis functions + int size ///< Kernel half-size ); Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 14669) +++ trunk/psModules/src/psmodules.h (revision 14671) @@ -79,4 +79,5 @@ #include #include +#include #include