Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 26156) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 26893) @@ -10,7 +10,9 @@ #include "pmSubtraction.h" #include "pmSubtractionKernels.h" +#include "pmSubtractionHermitian.h" +#include "pmSubtractionDeconvolve.h" +#include "pmSubtractionVisual.h" #define RINGS_BUFFER 10 // Buffer size for RINGS data - // Free function for pmSubtractionKernels @@ -27,4 +29,17 @@ psFree(kernels->solution1); psFree(kernels->solution2); + psFree(kernels->sampleStamps); +} + +// Free function for pmSubtractionPreCalcKernel +static void pmSubtractionKernelPreCalcFree(pmSubtractionKernelPreCalc *kernel) +{ + psFree(kernel->xKernel); + psFree(kernel->yKernel); + psFree(kernel->kernel); + + psFree(kernel->uCoords); + psFree(kernel->vCoords); + psFree(kernel->poly); } @@ -45,5 +60,5 @@ // Generate 1D convolution kernel for ISIS -static psVector *subtractionKernelISIS(float sigma, // Gaussian width +psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width int order, // Polynomial order int size // Kernel half-size @@ -57,4 +72,50 @@ 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) +psVector *pmSubtractionKernelHERM(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); + } + + return kernel; +} + +// Generate 1D convolution kernel for HERM (normalized for 2D) +psKernel *pmSubtractionKernelHERM_RADIAL(float sigma, // Gaussian width + int order, // Polynomial order + int size // Kernel half-size + ) +{ + psKernel *kernel = psKernelAlloc(-size, size, -size, size); // 2D Kernel + + // 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)); + + // generate 2D radial hermitian + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + float r = hypot(u, v) / sigma; + float z = -0.25*r*r; + kernel->kernel[v][u] = norm * p_pmSubtractionHermitianPolynomial(r, order) * exp(z); + } } @@ -96,5 +157,5 @@ kernels->preCalc->data[index] = NULL; kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); - + psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty"); psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); } @@ -110,15 +171,127 @@ } +bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, + int index, int size, int uOrder, int vOrder, float fwhm, + bool AlardLuptonStyle, bool forceZeroNull) +{ + // we have 4 cases here: + // 1) for odd functions, normalize the kernel by the maximum swing / Npix + // 2) for even functions, normalize the kernel to unity + // 3) for alard-lupton style normalization, subtract 1 from the 0,0 pixel for all even functions + // 4) for deconvolved hermitians, subtract 1 from the 0,0 pixel for the 0,0 function(s) + + // Calculate moments + double penalty = 0.0; // Moment, for penalty + double sum = 0.0, sum2 = 0.0; // Sum of kernel component + float min = INFINITY, max = -INFINITY; // Minimum and maximum kernel value + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + double value = preCalc->kernel->kernel[v][u]; + double value2 = PS_SQR(value); + sum += value; + sum2 += value2; + penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v))); + min = PS_MIN(value, min); + max = PS_MAX(value, max); + } + } + +#if 0 + fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty); +#endif + + bool zeroNull = false; // Zero out using the null position? + float scale2D = NAN; // Scaling for 2-D kernels + + if (AlardLuptonStyle) { + if (uOrder % 2 == 0 && vOrder % 2 == 0) { + // Even functions: normalise to unit sum and subtract null pixel so that sum is zero + scale2D = 1.0 / fabs(sum); + zeroNull = true; + } else { + // Odd functions: choose normalisation so that parameters have about the same strength as for even + // functions, no subtraction of null pixel because the sum is already (near) zero + scale2D = 1.0 / sqrt(sum2); + zeroNull = false; + } + } + + if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) { + zeroNull = true; + } + if (forceZeroNull) { + // Force rescaling and subtraction of null pixel even though the order doesn't indicate it's even + scale2D = 1.0 / fabs(sum); + zeroNull = true; + } + if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) { + // Odd function + scale2D = 1.0 / sqrt(sum2); + } + + float scale1D = sqrtf(scale2D); // Scaling for 1-D kernels + if (preCalc->xKernel) { + psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32)); + } + if (preCalc->yKernel) { + psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32)); + } + + psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32)); + penalty *= 1.0 / sum2; + + if (zeroNull) { + preCalc->kernel->kernel[0][0] -= 1.0; + } + +#if 0 + { + double sum = 0.0; // Sum of kernel component + float min = INFINITY, max = -INFINITY; // Minimum and maximum kernel value + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += preCalc->kernel->kernel[v][u]; + min = PS_MIN(preCalc->kernel->kernel[v][u], min); + max = PS_MAX(preCalc->kernel->kernel[v][u], max); + } + } + fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D); + } +#endif + + kernels->widths->data.F32[index] = fwhm; + 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 * penalty; + psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty"); + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty); + + return true; +} + pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, - const psVector *fwhms, const psVector *orders, - float penalty, pmSubtractionMode mode) -{ - 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); + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, psRegion bounds, 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 @@ -133,5 +306,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); @@ -141,5 +314,4 @@ // 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 @@ -147,60 +319,225 @@ 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] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel - psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(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))); - } + + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); + // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false); + } + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder, + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, psRegion bounds, 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; + num += (11 - gaussOrder - 1); // include all higher order radial terms + } + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, + spatialOrder, penalty, bounds, mode); // Kernels + psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS_RADIAL kernel: %s,%d --> %d elements", params, spatialOrder, num); + psFree(params); + + // 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 + // 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++) { + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); + } + } + for (int order = orders->data.S32[i] + 1; order < 11; order ++, index ++) { + // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian) + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true); + } + } + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder, + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, psRegion bounds, 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, bounds, mode); // 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 + 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++) { + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); + } + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder, + const psVector *fwhmsIN, const psVector *ordersIN, + float penalty, psRegion bounds, 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 += PS_SQR(gaussOrder + 1); + } + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, + spatialOrder, penalty, bounds, mode); // Kernels + psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num); + psFree(params); + + // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest) + // generate the Gaussian deconvolution kernel + # define DECONV_SIGMA 1.6 + psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECONV_SIGMA); + +# if 1 + psArray *deconKernels = psArrayAllocEmpty(100); +# endif + + // 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 + // Iterate over (u,v) order + for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { + for (int vOrder = 0; vOrder <= orders->data.S32[i]; vOrder++, index++) { + + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values + + // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel + psKernel *kernelTarget = preCalc->kernel; + preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel + + // XXX do we use Alard-Lupton normalization (last param true) or not? + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); + + // XXXX test demo that deconvolved kernel is valid +# if 1 + psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); + psArrayAdd (deconKernels, 100, kernelConv); + psFree (kernelConv); + + if (!uOrder && !vOrder){ + pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); } - - // 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)); - kernel->kernel[0][0] -= 1.0; - moment *= PS_SQR(sum); +# endif + } + } + } + +# if 1 + psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32); + for (int i = 0; i < deconKernels->n; i++) { + for (int j = 0; j <= i; j++) { + psImage *t1 = deconKernels->data[i]; + psImage *t2 = deconKernels->data[j]; + + double sum = 0.0; + for (int iy = 0; iy < t1->numRows; iy++) { + for (int ix = 0; ix < t1->numCols; ix++) { + sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix]; } - - -#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)); } - } - } + dot->data.F32[j][i] = sum; + dot->data.F32[i][j] = sum; + } + } + pmSubtractionVisualShowSubtraction (dot, NULL, NULL); + psFree (dot); + psFree (deconKernels); +# endif return kernels; @@ -212,5 +549,5 @@ pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, - int size, int spatialOrder, float penalty, + int size, int spatialOrder, float penalty, psRegion bounds, pmSubtractionMode mode) { @@ -224,9 +561,13 @@ kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); + kernels->uStop = NULL; + kernels->vStop = NULL; + kernels->xMin = bounds.x0; + kernels->xMax = bounds.x1; + kernels->yMin = bounds.y0; + kernels->yMax = bounds.y1; kernels->preCalc = psArrayAlloc(numBasisFunctions); kernels->penalty = penalty; kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); - kernels->uStop = NULL; - kernels->vStop = NULL; kernels->size = size; kernels->inner = 0; @@ -234,13 +575,77 @@ kernels->bgOrder = 0; kernels->mode = mode; - kernels->numCols = 0; - kernels->numRows = 0; kernels->solution1 = NULL; kernels->solution2 = NULL; + kernels->mean = NAN; + kernels->rms = NAN; + kernels->numStamps = 0; + kernels->sampleStamps = NULL; + + kernels->fSigResMean = NAN; + kernels->fSigResStdev = NAN; + kernels->fMaxResMean = NAN; + kernels->fMaxResStdev = NAN; + kernels->fMinResMean = NAN; + kernels->fMinResStdev = NAN; return kernels; } -pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, +pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(pmSubtractionKernelsType type, int uOrder, int vOrder, int size, float sigma) { + + pmSubtractionKernelPreCalc *preCalc = psAlloc(sizeof(pmSubtractionKernelPreCalc)); // Kernels, to return + psMemSetDeallocator(preCalc, (psFreeFunc)pmSubtractionKernelPreCalcFree); + + // 1D kernel realizations: + switch (type) { + case PM_SUBTRACTION_KERNEL_ISIS: + preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size); + preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size); + preCalc->uCoords = NULL; + preCalc->vCoords = NULL; + preCalc->poly = NULL; + break; + case PM_SUBTRACTION_KERNEL_HERM: + preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); + preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); + preCalc->uCoords = NULL; + preCalc->vCoords = NULL; + preCalc->poly = NULL; + break; + case PM_SUBTRACTION_KERNEL_RINGS: + // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure + // we allocate these vectors here, but leave the kernel generation to the main function + preCalc->xKernel = NULL; + preCalc->yKernel = NULL; + preCalc->kernel = NULL; + preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords + preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords + preCalc->poly = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial + return preCalc; + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + preCalc->kernel = pmSubtractionKernelHERM_RADIAL(sigma, uOrder, size); + preCalc->xKernel = NULL; + preCalc->yKernel = NULL; + preCalc->uCoords = NULL; + preCalc->vCoords = NULL; + preCalc->poly = NULL; + return preCalc; + default: + psAbort("programming error: invalid type for PreCalc kernel"); + } + + preCalc->kernel = psKernelAlloc(-size, size, -size, size); // 2D Kernel + + // generate 2D kernel from 1D realizations + for (int v = -size, y = 0; v <= size; v++, y++) { + for (int u = -size, x = 0; u <= size; u++, x++) { + preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel + } + } + + return preCalc; +} + +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, psRegion bounds, pmSubtractionMode mode) { @@ -251,5 +656,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty); psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", @@ -266,8 +671,8 @@ pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, const psVector *fwhms, const psVector *orders, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, - penalty, mode); // Kernels + penalty, bounds, mode); // Kernels if (!kernels) { return NULL; @@ -278,5 +683,5 @@ pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -299,5 +704,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder, @@ -370,5 +775,5 @@ pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty, - pmSubtractionMode mode) + psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -397,5 +802,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty); @@ -463,8 +868,8 @@ } -// Grid United with Normal Kernel +// Grid United with Normal Kernel [description: GUNK=ISIS(...)+POIS(...)] pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, const psVector *orders, int inner, float penalty, - pmSubtractionMode mode) + psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -479,5 +884,5 @@ pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, - penalty, mode); // Kernels + penalty, bounds, mode); // Kernels kernels->type = PM_SUBTRACTION_KERNEL_GUNK; psStringPrepend(&kernels->description, "GUNK="); @@ -495,5 +900,5 @@ // RINGS --- just what it says pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, - float penalty, pmSubtractionMode mode) + float penalty, psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -526,5 +931,5 @@ pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, - spatialOrder, penalty, mode); // The kernels + spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder, @@ -566,15 +971,15 @@ for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) { - psArray *data = psArrayAlloc(3); // Container for data - psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords - psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords - psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0); double moment = 0.0; // Moment, for penalty if (i == 0) { // Central pixel is easy - uCoords->data.S32[0] = vCoords->data.S32[0] = 0; - poly->data.F32[0] = 1.0; - uCoords->n = vCoords->n = poly->n = 1; + preCalc->uCoords->data.S32[0] = 0; + preCalc->vCoords->data.S32[0] = 0; + preCalc->poly->data.F32[0] = 1.0; + preCalc->uCoords->n = 1; + preCalc->vCoords->n = 1; + preCalc->poly->n = 1; radiusLast = 0; moment = 0.0; @@ -594,15 +999,15 @@ float polyVal = uPoly * vPoly; // Value of polynomial if (polyVal != 0) { // No point adding it otherwise - uCoords->data.S32[j] = u; - vCoords->data.S32[j] = v; - poly->data.F32[j] = polyVal; + preCalc->uCoords->data.S32[j] = u; + preCalc->vCoords->data.S32[j] = v; + preCalc->poly->data.F32[j] = polyVal; norm += polyVal; - moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v)); - - psVectorExtend(uCoords, RINGS_BUFFER, 1); - psVectorExtend(vCoords, RINGS_BUFFER, 1); - psVectorExtend(poly, RINGS_BUFFER, 1); + moment += PS_SQR(polyVal) * PS_SQR(PS_SQR(u) + PS_SQR(v)); + + psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->poly, RINGS_BUFFER, 1); psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n", - u, v, poly->data.F32[j]); + u, v, preCalc->poly->data.F32[j]); j++; } @@ -612,25 +1017,28 @@ // Normalise kernel component to unit sum if (uOrder % 2 == 0 && vOrder % 2 == 0) { - psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); + psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); // Add subtraction of 0,0 component to preserve photometric scaling - uCoords->data.S32[j] = 0; - vCoords->data.S32[j] = 0; - poly->data.F32[j] = -1.0; - psVectorExtend(uCoords, RINGS_BUFFER, 1); - psVectorExtend(vCoords, RINGS_BUFFER, 1); - psVectorExtend(poly, RINGS_BUFFER, 1); + preCalc->uCoords->data.S32[j] = 0; + preCalc->vCoords->data.S32[j] = 0; + preCalc->poly->data.F32[j] = -1.0; + psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->poly, RINGS_BUFFER, 1); } else { norm = powf(size, uOrder) * powf(size, vOrder); - psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); + psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); } - moment /= norm; + moment /= PS_SQR(norm); } - psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", uCoords->n); - - kernels->preCalc->data[index] = data; + psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", preCalc->uCoords->n); + + kernels->preCalc->data[index] = preCalc; kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); + if (!isfinite(kernels->penalties->data.F32[index])) { + psAbort ("invalid penalty"); + } psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, @@ -645,20 +1053,26 @@ pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *fwhms, const psVector *orders, int inner, - int binning, int ringsOrder, float penalty, + int binning, int ringsOrder, float penalty, psRegion bounds, pmSubtractionMode mode) { switch (type) { case PM_SUBTRACTION_KERNEL_POIS: - return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode); + return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_ISIS: - return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode); + return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode); + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, bounds, mode); + case PM_SUBTRACTION_KERNEL_HERM: + return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode); + case PM_SUBTRACTION_KERNEL_DECONV_HERM: + return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_SPAM: - return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_FRIES: - return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode); + return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_GUNK: - return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode); + return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, bounds, mode); case PM_SUBTRACTION_KERNEL_RINGS: - return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode); + return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); @@ -696,5 +1110,5 @@ pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder, - pmSubtractionMode mode) + psRegion bounds, pmSubtractionMode mode) { PS_ASSERT_STRING_NON_EMPTY(description, NULL); @@ -715,36 +1129,40 @@ float penalty = 0.0; // Penalty for wideness - if (strncmp(description, "ISIS", 4) == 0) { - // XXX Support for GUNK - 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 = PM_SUBTRACTION_KERNEL_RINGS; - char *ptr = (char*)description + 6; + // currently known descriptions: + // ISIS(...), ISIS_RADIAL(...), HERM(...), DECONV_HERM(...), POIS(...), SPAM(...), + // FRIES(...), GUNK=ISIS(...)+POIS(...), RINGS(...), + // the descriptive name is the set of characters before the ( + + type = pmSubtractionKernelsTypeFromString (description); + char *ptr = strchr(description, '(') + 1; + psAssert (ptr, "description is missing kernel parameters"); + + switch (type) { + case PM_SUBTRACTION_KERNEL_ISIS: + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + case PM_SUBTRACTION_KERNEL_HERM: + case PM_SUBTRACTION_KERNEL_DECONV_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); + break; + case PM_SUBTRACTION_KERNEL_RINGS: PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt); @@ -752,32 +1170,56 @@ 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); -} - - + break; + default: + psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported."); + } + + return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, + ringsOrder, penalty, bounds, mode); +} + + +// the input string can either be just the name or the description string. Currently known +// descriptions: ISIS(...), ISIS_RADIAL(...), HERM(...), DECONV_HERM(...), POIS(...), +// SPAM(...), FRIES(...), GUNK=ISIS(...)+POIS(...), RINGS(...), pmSubtractionKernelsType pmSubtractionKernelsTypeFromString(const char *type) { - if (strcasecmp(type, "POIS") == 0) { + // for a bare name (ISIS, HERM), use the full string length. + // otherwise, use the length up to the first '(' + int nameLength = strlen(type); + char *ptr = strchr(type, '('); + if (ptr) { + nameLength = ptr - type; + } + + if (strncasecmp(type, "POIS", nameLength) == 0) { return PM_SUBTRACTION_KERNEL_POIS; } - if (strcasecmp(type, "ISIS") == 0) { + if (strncasecmp(type, "ISIS", nameLength) == 0) { return PM_SUBTRACTION_KERNEL_ISIS; } - if (strcasecmp(type, "SPAM") == 0) { + if (strncasecmp(type, "ISIS_RADIAL", nameLength) == 0) { + return PM_SUBTRACTION_KERNEL_ISIS_RADIAL; + } + if (strncasecmp(type, "HERM", nameLength) == 0) { + return PM_SUBTRACTION_KERNEL_HERM; + } + if (strncasecmp(type, "DECONV_HERM", nameLength) == 0) { + return PM_SUBTRACTION_KERNEL_DECONV_HERM; + } + if (strncasecmp(type, "SPAM", nameLength) == 0) { return PM_SUBTRACTION_KERNEL_SPAM; } - if (strcasecmp(type, "FRIES") == 0) { + if (strncasecmp(type, "FRIES", nameLength) == 0) { return PM_SUBTRACTION_KERNEL_FRIES; } - if (strcasecmp(type, "GUNK") == 0) { + if (strncasecmp(type, "GUNK", nameLength) == 0) { return PM_SUBTRACTION_KERNEL_GUNK; } - if (strcasecmp(type, "RINGS") == 0) { + // note that GUNK has a somewhat different description + if (strncasecmp(type, "GUNK=ISIS", nameLength) == 0) { + return PM_SUBTRACTION_KERNEL_GUNK; + } + if (strncasecmp(type, "RINGS", nameLength) == 0) { return PM_SUBTRACTION_KERNEL_RINGS; } @@ -810,8 +1252,11 @@ out->bgOrder = in->bgOrder; out->mode = in->mode; - out->numCols = in->numCols; - out->numRows = in->numRows; + out->xMin = in->xMin; + out->xMax = in->xMax; + out->yMin = in->yMin; + out->yMax = in->yMax; out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL; out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL; + out->sampleStamps = psMemIncrRefCounter(in->sampleStamps); return out;