Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14542) +++ 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