Index: /branches/pap/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 26716) +++ /branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 26717) @@ -188,5 +188,11 @@ if (normalise) { // Put in the normalisation component - kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels)); +// kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels)); + double value = wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels); + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + kernel->kernel[v][u] += value * kernels->norm->kernel[v][u]; + } + } } @@ -754,7 +760,9 @@ case PM_SUBTRACTION_MODE_1: stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); + stamp->normConv = p_pmSubtractionConvolveStampPrecalc(stamp->image1, kernels->norm); break; case PM_SUBTRACTION_MODE_2: stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint); + stamp->normConv = p_pmSubtractionConvolveStampPrecalc(stamp->image2, kernels->norm); break; case PM_SUBTRACTION_MODE_UNSURE: @@ -762,4 +770,5 @@ stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint); + stamp->normConv = p_pmSubtractionConvolveStampPrecalc(stamp->image1, kernels->norm); break; default: @@ -908,4 +917,7 @@ psFree(stamp->vector); stamp->vector = NULL; + + psFree(stamp->normConv); + stamp->normConv = NULL; } else { numGood++; Index: /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 26716) +++ /branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 26717) @@ -229,4 +229,5 @@ const psArray *convolutions1, // Convolutions of image 1 for each kernel const psArray *convolutions2, // Convolutions of image 2 for each kernel + const psKernel *normConv, // Normalisation, convolved const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values @@ -370,8 +371,13 @@ double sumB = 0.0; // Sum of B products (for matrix, background) double sumI2 = 0.0; // Sum of I_2 (for vector, background) + + double sumAN = 0.0; // Matrix: A*norm + double sumBN = 0.0; // Matrix: B*norm + for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { double a = iConv1->kernel[y][x]; double b = iConv2->kernel[y][x]; + double n = normConv->kernel[y][x]; float i1 = image1->kernel[y][x]; float i2 = image2->kernel[y][x]; @@ -381,4 +387,7 @@ double ai1 = a * i1; double bi1 = b * i1; + + double an = a * n; + double bn = b * n; if (weight) { @@ -391,4 +400,7 @@ b *= wtVal; i2 *= wtVal; + + an *= wtVal; + bn *= wtVal; } if (window) { @@ -401,4 +413,7 @@ b *= wtVal; i2 *= wtVal; + + an *= wtVal; + bn *= wtVal; } sumAI2 += ai2; @@ -409,4 +424,7 @@ sumB += b; sumI2 += i2; + + sumAN += an; + sumBN += bn; } } @@ -420,9 +438,16 @@ double b = sumB * poly[iTerm]; + double an = sumAN * poly[iTerm]; + double bn = sumBN * poly[iTerm]; + if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = ai1; - matrix->data.F64[normIndex][iIndex] = ai1; - matrix->data.F64[iIndex + numParams][normIndex] = bi1; - matrix->data.F64[normIndex][iIndex + numParams] = bi1; +// matrix->data.F64[iIndex][normIndex] = ai1; +// matrix->data.F64[normIndex][iIndex] = ai1; +// matrix->data.F64[iIndex + numParams][normIndex] = bi1; +// matrix->data.F64[normIndex][iIndex + numParams] = bi1; + matrix->data.F64[iIndex][normIndex] = an; + matrix->data.F64[normIndex][iIndex] = an; + matrix->data.F64[iIndex + numParams][normIndex] = bn; + matrix->data.F64[normIndex][iIndex + numParams] = bn; } if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { @@ -453,4 +478,9 @@ double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window + + double sumN = 0.0; // Matrix: bg*norm + double sumN2 = 0.0; // Matrix: norm*norm + double sumNI2 = 0.0; // Vector: I2*norm + for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -461,4 +491,8 @@ double one = 1.0; double i1i2 = i1 * i2; + + double n = normConv->kernel[y][x]; + double n2 = n * n; + double ni2 = n * i2; if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { @@ -474,4 +508,8 @@ i2 *= wtVal; i1i2 *= wtVal; + + n *= wtVal; + n2 *= wtVal; + ni2 *= wtVal; } if (window) { @@ -482,4 +520,8 @@ i2 *= wtVal; i1i2 *= wtVal; + + n *= wtVal; + n2 *= wtVal; + ni2 *= wtVal; } sumI1 += i1; @@ -488,4 +530,8 @@ sumI2 += i2; sumI1I2 += i1i2; + + sumN += n; + sumN2 += n2; + sumNI2 += ni2; } } @@ -496,6 +542,9 @@ if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumI1I1; - vector->data.F64[normIndex] = sumI1I2; +// matrix->data.F64[normIndex][normIndex] = sumI1I1; +// vector->data.F64[normIndex] = sumI1I2; + matrix->data.F64[normIndex][normIndex] = sumN2; + vector->data.F64[normIndex] = sumNI2; + } if (mode & PM_SUBTRACTION_EQUATION_BG) { @@ -504,6 +553,8 @@ } if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[bgIndex][normIndex] = sumI1; - matrix->data.F64[normIndex][bgIndex] = sumI1; +// matrix->data.F64[bgIndex][normIndex] = sumI1; +// matrix->data.F64[normIndex][bgIndex] = sumI1; + matrix->data.F64[bgIndex][normIndex] = sumN; + matrix->data.F64[normIndex][bgIndex] = sumN; } return true; @@ -729,5 +780,6 @@ stamp->image1, stamp->image2, weight, window, stamp->convolutions1, stamp->convolutions2, - kernels, polyValues, footprint, stamps->normWindow, mode); + stamp->normConv, kernels, polyValues, footprint, + stamps->normWindow, mode); break; default: Index: /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c (revision 26716) +++ /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c (revision 26717) @@ -15,4 +15,6 @@ #define RINGS_BUFFER 10 // Buffer size for RINGS data + +#define NORM_WIDTH 7.5 // Free function for pmSubtractionKernels @@ -30,4 +32,5 @@ psFree(kernels->solution2); psFree(kernels->sampleStamps); + psFree(kernels->norm); } @@ -242,5 +245,10 @@ if (zeroNull) { - preCalc->kernel->kernel[0][0] -= 1.0; + // preCalc->kernel->kernel[0][0] -= 1.0; + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + preCalc->kernel->kernel[v][u] -= kernels->norm->kernel[v][u]; + } + } } @@ -582,4 +590,19 @@ kernels->fMinResMean = NAN; kernels->fMinResStdev = NAN; + + kernels->norm = psKernelAlloc(-size, size, -size, size); + double sumNorm = 0.0; + float sigma = NORM_WIDTH / 2.35; + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + float z = 0.5 * (PS_SQR(x) + PS_SQR(y)) / PS_SQR(sigma); + sumNorm += kernels->norm->kernel[y][x] = expf(-z); + } + } + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + kernels->norm->kernel[y][x] /= sumNorm; + } + } return kernels; Index: /branches/pap/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionKernels.h (revision 26716) +++ /branches/pap/psModules/src/imcombine/pmSubtractionKernels.h (revision 26717) @@ -12,5 +12,5 @@ PM_SUBTRACTION_KERNEL_ISIS_RADIAL, ///< ISIS + higher-order radial Hermitians PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels - PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved Hermitian polynomial kernels + PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved 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 @@ -37,4 +37,5 @@ psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM) + psKernel *norm; ///< Normalisation component float penalty; ///< Penalty for wideness psVector *penalties; ///< Penalty for each kernel component @@ -49,22 +50,22 @@ float mean, rms; ///< Mean and RMS of chi^2 from stamps int numStamps; ///< Number of good stamps - float fSigResMean; ///< mean fractional stdev of residuals - float fSigResStdev; ///< stdev of fractional stdev of residuals - float fMaxResMean; ///< mean fractional positive swing in residuals - float fMaxResStdev; ///< stdev of fractional positive swing in residuals - float fMinResMean; ///< mean fractional negative swing in residuals - float fMinResStdev; ///< stdev of fractional negative swing in residuals - psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations + float fSigResMean; ///< mean fractional stdev of residuals + float fSigResStdev; ///< stdev of fractional stdev of residuals + float fMaxResMean; ///< mean fractional positive swing in residuals + float fMaxResStdev; ///< stdev of fractional positive swing in residuals + float fMinResMean; ///< mean fractional negative swing in residuals + float fMinResStdev; ///< stdev of fractional negative swing in residuals + psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations } pmSubtractionKernels; // pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures typedef struct { - psVector *uCoords; // used by RINGS - psVector *vCoords; // used by RINGS - psVector *poly; // used by RINGS - - psVector *xKernel; // used by ISIS, HERM, DECONV_HERM - psVector *yKernel; // used by ISIS, HERM, DECONV_HERM - psKernel *kernel; // used by ISIS, HERM, DECONV_HERM + psVector *uCoords; // used by RINGS + psVector *vCoords; // used by RINGS + psVector *poly; // used by RINGS + + psVector *xKernel; // used by ISIS, HERM, DECONV_HERM + psVector *yKernel; // used by ISIS, HERM, DECONV_HERM + psKernel *kernel; // used by ISIS, HERM, DECONV_HERM } pmSubtractionKernelPreCalc; @@ -168,8 +169,8 @@ pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc( pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated) - int uOrder, ///< order in x-direction - int vOrder, ///< order in x-direction - int size, ///< Half-size of the kernel - float sigma ///< sigma of gaussian kernel + int uOrder, ///< order in x-direction + int vOrder, ///< order in x-direction + int size, ///< Half-size of the kernel + float sigma ///< sigma of gaussian kernel ); @@ -202,9 +203,9 @@ /// Generate ISIS + RADIAL_HERM kernels pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(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 - float penalty, ///< Penalty for wideness - pmSubtractionMode mode ///< Mode for subtraction + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< Polynomial order of gaussians + float penalty, ///< Penalty for wideness + pmSubtractionMode mode ///< Mode for subtraction ); @@ -220,9 +221,9 @@ /// Generate DECONV_HERM kernels pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(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 + 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 ); Index: /branches/pap/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionStamps.c (revision 26716) +++ /branches/pap/psModules/src/imcombine/pmSubtractionStamps.c (revision 26717) @@ -312,4 +312,5 @@ stamp->vector = NULL; stamp->norm = NAN; + stamp->normConv = NULL; return stamp; Index: /branches/pap/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/pap/psModules/src/imcombine/pmSubtractionStamps.h (revision 26716) +++ /branches/pap/psModules/src/imcombine/pmSubtractionStamps.h (revision 26717) @@ -83,4 +83,5 @@ psVector *vector; ///< Least-squares vector, or NULL double norm; ///< Normalisation difference + psKernel *normConv; ///< Convolution of normalisation pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp;