Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26472) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26473) @@ -10,7 +10,7 @@ #include "pmSubtraction.h" #include "pmSubtractionKernels.h" +#include "pmSubtractionHermitian.h" #define RINGS_BUFFER 10 // Buffer size for RINGS data - // Free function for pmSubtractionKernels @@ -57,4 +57,26 @@ for (int i = 0, x = -size; x <= size; i++, x++) { kernel->data.F32[i] = norm * power(x, order) * expf(expNorm * PS_SQR(x)); + } + + return kernel; +} + +// Generate 1D convolution kernel for HERM (normalized for 2D) +static psVector *subtractionKernelHERM(float sigma, // Gaussian width + int order, // Polynomial order + int size // Kernel half-size + ) +{ + int fullSize = 2 * size + 1; // Full size of kernel + psVector *kernel = psVectorAlloc(fullSize, PS_TYPE_F32); // Kernel to return + + // for now, we are only allowing equal orders and sigmas in X and Y + float nf = exp(lgamma(order + 1)); + float norm = 1.0 / sqrt(nf*sigma*sqrt(M_2_PI)); + + for (int i = 0, x = -size; x <= size; i++, x++) { + float xf = x / sigma; + float z = -0.25*xf*xf; + kernel->data.F32[i] = norm * p_pmSubtractionHermitianPolynomial(xf, order) * exp(z); } @@ -189,4 +211,113 @@ #if 0 + double sum = 0.0; // Sum of kernel component + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += kernel->kernel[v][u]; + } + } + fprintf(stderr, "%d sum: %lf\n", index, sum); +#endif + + kernels->widths->data.F32[index] = fwhms->data.F32[i]; + 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; + kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); + + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, + fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); + } + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder, + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, pmSubtractionMode mode) +{ + PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); + PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL); + PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL); + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + + // check the requested fwhm values: any values <= 0.0 should be dropped + psVector *fwhms = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32); + psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32); + for (int i = 0; i < fwhmsIN->n; i++) { + if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; + psVectorAppend(fwhms, fwhmsIN->data.F32[i]); + psVectorAppend(orders, ordersIN->data.S32[i]); + } + + int numGaussians = fwhms->n; // Number of Gaussians + + int num = 0; // Number of basis functions + psString params = NULL; // List of parameters + for (int i = 0; i < numGaussians; i++) { + int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian + psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); + num += (gaussOrder + 1) * (gaussOrder + 2) / 2; + } + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", params, spatialOrder, num); + psFree(params); + + // Set the kernel parameters + int fullSize = 2 * size + 1; // Full size of kernels + 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 + // 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++) { + psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values + psVector *xKernel = preCalc->data[0] = subtractionKernelHERM(sigma, uOrder, size); // x Kernel + psVector *yKernel = preCalc->data[1] = subtractionKernelHERM(sigma, vOrder, size); // y Kernel + psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size); // Kernel + + // Calculate moments + double moment = 0.0; // Moment, for penalty + for (int v = -size, y = 0; v <= size; v++, y++) { + for (int u = -size, x = 0; u <= size; u++, x++) { + double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel + kernel->kernel[v][u] = value; + moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); + } + } + + // Normalise sum of kernel component to unity for even functions + if (uOrder % 2 == 0 && vOrder % 2 == 0) { + double sum = 0.0; // Sum of kernel component + for (int v = 0; v < fullSize; v++) { + for (int u = 0; u < fullSize; u++) { + sum += xKernel->data.F32[u] * yKernel->data.F32[v]; + } + } + sum = 1.0 / sqrt(sum); + psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); + psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); + psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); + +#if 1 + fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum,kernel->kernel[0][0]); +#endif + + kernel->kernel[0][0] -= 1.0; + moment *= PS_SQR(sum); + } + + +#if 1 double sum = 0.0; // Sum of kernel component for (int v = -size; v <= size; v++) { @@ -662,4 +793,6 @@ case PM_SUBTRACTION_KERNEL_ISIS: return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode); + case PM_SUBTRACTION_KERNEL_HERM: + return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode); case PM_SUBTRACTION_KERNEL_SPAM: return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); @@ -724,34 +857,41 @@ float penalty = 0.0; // Penalty for wideness - if (strncmp(description, "ISIS", 4) == 0) { - // XXX Support for GUNK + // ISIS and HERM have the same description layout + if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM")) { + // XXX Support for GUNK (not yet supported) if (strstr(description, "+POIS")) { type = PM_SUBTRACTION_KERNEL_GUNK; psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description); - } else { - type = PM_SUBTRACTION_KERNEL_ISIS; - char *ptr = (char*)description + 5; // Eat "ISIS(" - PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); - - // Count the number of Gaussians - int numGauss = 0; - for (char *string = ptr; string; string = strchr(string + 1, '(')) { - numGauss++; - } - - fwhms = psVectorAlloc(numGauss, PS_TYPE_F32); - orders = psVectorAlloc(numGauss, PS_TYPE_S32); - - for (int i = 0; i < numGauss; i++) { - ptr++; // Eat the '(' - PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234," - PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)" - } - - ptr++; // Eat ',' - PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); - penalty = parseStringFloat(ptr); - } - } else if (strncmp(description, "RINGS", 5) == 0) { + } + + type = pmSubtractionKernelsTypeFromString (description); + psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must be ISIS or HERM"); + + char *ptr = (char*)description + 5; // Eat "ISIS(" or "HERM(" + PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); + + // Count the number of Gaussians + int numGauss = 0; + for (char *string = ptr; string; string = strchr(string + 1, '(')) { + numGauss++; + } + + fwhms = psVectorAlloc(numGauss, PS_TYPE_F32); + orders = psVectorAlloc(numGauss, PS_TYPE_S32); + + for (int i = 0; i < numGauss; i++) { + ptr++; // Eat the '(' + PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234," + PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)" + } + + ptr++; // Eat ',' + PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); + penalty = parseStringFloat(ptr); + + return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); + } + + if (strncmp(description, "RINGS", 5) == 0) { type = PM_SUBTRACTION_KERNEL_RINGS; char *ptr = (char*)description + 6; @@ -761,8 +901,8 @@ PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); - } else { - psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported."); - } - + return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); + } + + psAbort("Deciphering kernels other than ISIS, HERM and RINGS is not currently supported."); return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, @@ -778,4 +918,7 @@ if (strcasecmp(type, "ISIS") == 0) { return PM_SUBTRACTION_KERNEL_ISIS; + } + if (strcasecmp(type, "HERM") == 0) { + return PM_SUBTRACTION_KERNEL_HERM; } if (strcasecmp(type, "SPAM") == 0) { Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26472) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26473) @@ -10,4 +10,5 @@ PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction @@ -30,8 +31,8 @@ psString description; ///< Description of the kernel parameters long num; ///< Number of kernel components (not including the spatial ones) - psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) - psVector *widths; ///< Gaussian FWHMs (ISIS) + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS or HERM) + psVector *widths; ///< Gaussian FWHMs (ISIS or HERM) psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS or HERM) float penalty; ///< Penalty for wideness psVector *penalties; ///< Penalty for each kernel component @@ -64,4 +65,9 @@ PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \ } \ + if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_HERM) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \ + } \ if ((KERNELS)->uStop || (KERNELS)->vStop) { \ PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \ @@ -142,4 +148,13 @@ ); +/// Generate HERM kernels +pmSubtractionKernels *pmSubtractionKernelsHERM(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< order of hermitian polynomials + float penalty, ///< Penalty for wideness + pmSubtractionMode mode ///< Mode for subtraction + ); + /// Generate SPAM kernels pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel