Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c (revision 29164) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c (revision 29165) @@ -27,7 +27,7 @@ // Calculate the least-squares matrix and vector -static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated - psVector *vector, // Least-squares vector, updated - double *norm, // Normalisation, updated +static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated + double normValue, // Normalisation, supplied const psKernel *input, // Input image (target) const psKernel *reference, // Reference image (convolution source) @@ -37,8 +37,5 @@ const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - int normWindow1, // Window (half-)size for normalisation measurement - int normWindow2, // Window (half-)size for normalisation measurement - const pmSubtractionEquationCalculationMode mode + int footprint // (Half-)Size of stamp ) { @@ -50,13 +47,21 @@ int numKernels = kernels->num; // Number of kernels - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background int spatialOrder = kernels->spatialOrder; // Order of spatial variation int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms double poly[numPoly]; // Polynomial terms double poly2[numPoly][numPoly]; // Polynomial-polynomial values + int numParams = numKernels * numPoly; + + psAssert(matrix && + matrix->type.type == PS_TYPE_F64 && + matrix->numCols == numParams && + matrix->numRows == numParams, + "Least-squares matrix is bad."); + psAssert(vector && + vector->type.type == PS_TYPE_F64 && + vector->n == numParams, + "Least-squares vector is bad."); // Evaluate polynomial-polynomial terms - // XXX we can skip this if we are not calculating kernel coeffs for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { @@ -84,6 +89,4 @@ // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] - // normalization - // bg 0, bg 1, bg 2 (only 0 is currently used?) for (int i = 0; i < numKernels; i++) { @@ -107,18 +110,15 @@ // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double value = sumCC * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex] = value; - matrix->data.F64[jIndex][iIndex] = value; - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double value = sumCC * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex] = value; + matrix->data.F64[jIndex][iIndex] = value; + } + } } double sumRC = 0.0; // Sum of the reference-convolution products double sumIC = 0.0; // Sum of the input-convolution products - double sumC = 0.0; // Sum of the convolution for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -128,10 +128,8 @@ double ic = in * conv; double rc = ref * conv; - double c = conv; if (weight) { float wtVal = weight->kernel[y][x]; ic *= wtVal; rc *= wtVal; - c *= wtVal; } if (window) { @@ -139,97 +137,14 @@ ic *= winVal; rc *= winVal; - c *= winVal; } sumIC += ic; sumRC += rc; - sumC += c; - } - } + } + } + // Spatial variation for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - double normTerm = sumRC * poly[iTerm]; - double bgTerm = sumC * poly[iTerm]; - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = normTerm; - matrix->data.F64[normIndex][iIndex] = normTerm; - } - if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = bgTerm; - matrix->data.F64[bgIndex][iIndex] = bgTerm; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = sumIC * poly[iTerm]; - if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { - // subtract norm * sumRC * poly[iTerm] - psAssert (kernels->solution1, "programming error: define solution first!"); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation - vector->data.F64[iIndex] -= norm * normTerm; - } - } - } - } - - double sumRR = 0.0; // Sum of the reference product - double sumIR = 0.0; // Sum of the input-reference product - double sum1 = 0.0; // Sum of the background - double sumR = 0.0; // Sum of the reference - double sumI = 0.0; // Sum of the input - double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double in = input->kernel[y][x]; - double ref = reference->kernel[y][x]; - double ir = in * ref; - double rr = PS_SQR(ref); - double one = 1.0; - - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { - normI1 += ref; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { - normI2 += in; - } - - if (weight) { - float wtVal = weight->kernel[y][x]; - rr *= wtVal; - ir *= wtVal; - in *= wtVal; - ref *= wtVal; - one *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - rr *= winVal; - ir *= winVal; - in *= winVal; - ref *= winVal; - one *= winVal; - } - sumRR += rr; - sumIR += ir; - sumR += ref; - sumI += in; - sum1 += one; - } - } - - *norm = normI2 / normI1; - - fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); - - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumRR; - vector->data.F64[normIndex] = sumIR; - // subtract sum over kernels * kernel solution - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI; - } - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[normIndex][bgIndex] = sumR; - matrix->data.F64[bgIndex][normIndex] = sumR; + vector->data.F64[iIndex] = (sumIC - normValue*sumRC) * poly[iTerm]; + } } @@ -255,7 +170,7 @@ // Calculate the least-squares matrix and vector for dual convolution -static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated - psVector *vector, // Least-squares vector, updated - double *norm, // Normalisation, updated +static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated + double normValue, // Normalisation, updated const psKernel *image1, // Image 1 const psKernel *image2, // Image 2 @@ -266,13 +181,8 @@ const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - int normWindow1, // Window (half-)size for normalisation measurement - int normWindow2, // Window (half-)size for normalisation measurement - const pmSubtractionEquationCalculationMode mode + int footprint // (Half-)Size of stamp ) { int numKernels = kernels->num; // Number of kernels - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background int spatialOrder = kernels->spatialOrder; // Order of spatial variation int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms @@ -280,8 +190,7 @@ double poly2[numPoly][numPoly]; // Polynomial-polynomial values - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters - int numParams2 = numKernels * numPoly; // Number of additional parameters for dual - int numDual = numParams + numParams2; // Total number of parameters for dual + int numParams = numKernels * numPoly; // Number of regular parameters + int numParams2 = numKernels * numPoly; // Number of additional parameters for dual + int numDual = numParams + numParams2; // Total number of parameters for dual psAssert(matrix && @@ -317,4 +226,9 @@ matrix->data.F64[i][i] = 1.0; } + + // the order of the elements in the matrix and vector is: + // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0] + // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] + // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] for (int i = 0; i < numKernels; i++) { @@ -352,23 +266,23 @@ // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double aa = sumAA * poly2[iTerm][jTerm]; - double bb = sumBB * poly2[iTerm][jTerm]; - double ab = sumAB * poly2[iTerm][jTerm]; - - matrix->data.F64[iIndex][jIndex] = aa; - matrix->data.F64[jIndex][iIndex] = aa; - - matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; - matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; - - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double aa = sumAA * poly2[iTerm][jTerm]; + double bb = sumBB * poly2[iTerm][jTerm]; + double ab = sumAB * poly2[iTerm][jTerm]; + + matrix->data.F64[iIndex][jIndex] = aa; + matrix->data.F64[jIndex][iIndex] = aa; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; + matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; + + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } + + // we need to calculate the lower-diagonal AB elements since they are not symmetric for A <-> B for (int j = 0; j < i; j++) { psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j @@ -388,12 +302,10 @@ // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double ab = sumAB * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double ab = sumAB * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } } } @@ -402,8 +314,5 @@ double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) - double sumA = 0.0; // Sum of A (for matrix, background) double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) - double sumB = 0.0; // Sum of B products (for matrix, background) - double sumI2 = 0.0; // Sum of I_2 (for vector, background) for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -424,7 +333,4 @@ ai1 *= wtVal; bi1 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; } if (window) { @@ -434,15 +340,9 @@ ai1 *= wtVal; bi1 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; } sumAI2 += ai2; sumBI2 += bi2; sumAI1 += ai1; - sumA += a; sumBI1 += bi1; - sumB += b; - sumI2 += i2; } } @@ -452,95 +352,8 @@ double bi2 = sumBI2 * poly[iTerm]; double ai1 = sumAI1 * poly[iTerm]; - double a = sumA * poly[iTerm]; double bi1 = sumBI1 * poly[iTerm]; - double b = sumB * 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; - } - if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = a; - matrix->data.F64[bgIndex][iIndex] = a; - matrix->data.F64[iIndex + numParams][bgIndex] = b; - matrix->data.F64[bgIndex][iIndex + numParams] = b; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = ai2; - vector->data.F64[iIndex + numParams] = bi2; - if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { - // subtract norm * sumRC * poly[iTerm] - psAssert (kernels->solution1, "programming error: define solution first!"); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation - vector->data.F64[iIndex] -= norm * ai1; - vector->data.F64[iIndex + numParams] -= norm * bi1; - } - } - } - } - - double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) - double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) - double sum1 = 0.0; // Sum of 1 (for matrix, background-background) - double sumI2 = 0.0; // Sum of I_2 (for vector, background) - 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 - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double i1 = image1->kernel[y][x]; - double i2 = image2->kernel[y][x]; - - double i1i1 = i1 * i1; - double one = 1.0; - double i1i2 = i1 * i2; - - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { - normI1 += i1; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { - normI2 += i2; - } - - if (weight) { - float wtVal = weight->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } - sumI1 += i1; - sumI1I1 += i1i1; - sum1 += one; - sumI2 += i2; - sumI1I2 += i1i2; - } - } - - *norm = normI2 / normI1; - fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); - - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumI1I1; - vector->data.F64[normIndex] = sumI1I2; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI2; - } - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[bgIndex][normIndex] = sumI1; - matrix->data.F64[normIndex][bgIndex] = sumI1; + vector->data.F64[iIndex] = ai2 - normValue * ai1; + vector->data.F64[iIndex + numParams] = bi2 - normValue * bi1; + } } @@ -565,12 +378,14 @@ } -#if 1 // Add in penalty term to least-squares vector bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term psVector *vector, // Vector to which to add in penalty term const pmSubtractionKernels *kernels, // Kernel parameters - float norm // Normalisation + float normSquare1, // Normalisation for image 1 + float normSquare2 // Normalisation for image 2 ) { + psAssert (kernels->mode == PM_SUBTRACTION_MODE_DUAL, "only use penalties for dual convolution"); + if (kernels->penalty == 0.0) { return true; @@ -589,6 +404,5 @@ // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] - // [norm] - // [bg] + // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] @@ -600,21 +414,19 @@ // Contribution to chi^2: a_i^2 P_i psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty"); - fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]); - matrix->data.F64[index][index] += norm * penalties1->data.F32[i]; - if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], - matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]); - matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i]; - // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) - // penalties scale with second moments - // - } - } - } - } - - return true; -} -# endif + fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]); + matrix->data.F64[index][index] += normSquare1 * penalties1->data.F32[i]; + + fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], + matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]); + matrix->data.F64[index + numParams][index + numParams] += normSquare2 * penalties2->data.F32[i]; + // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) + // penalties scale with second moments + // + } + } + } + + return true; +} ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -647,12 +459,5 @@ int index, bool wantDual) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); - PS_ASSERT_INT_POSITIVE(index, NAN); -#endif - psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, @@ -662,11 +467,5 @@ double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); -#endif - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation return kernels->solution1->data.F64[normIndex]; @@ -676,11 +475,5 @@ const psImage *polyValues) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); -#endif - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); @@ -690,4 +483,88 @@ // Public functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmSubtractionCalculateNormalization( + pmSubtractionStampList *stamps, + const pmSubtractionMode mode) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + psTimerStart("pmSubtractionCalculateNormalization"); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + + int footprint = stamps->footprint; // Half-size of stamps + + // Loop over each stamp and calculate its normalization factor + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; + if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; + + // XXX skip this if we have already calculated it? (stamp->norm does not change, just the median statistic) + // XXX maybe not: the star may have changed for a given stamp -- only if norm is reset to NAN can we do this + if (mode == PM_SUBTRACTION_MODE_2) { + pmSubtractionCalculateNormalizationStamp(stamp, stamp->image2, stamp->image1, footprint, stamps->normWindow2, stamps->normWindow1); + } else { + pmSubtractionCalculateNormalizationStamp(stamp, stamp->image1, stamp->image2, footprint, stamps->normWindow1, stamps->normWindow2); + } + psVectorAppend(norms, stamp->norm); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(norms); + return false; + } + stamps->normValue = stats->robustMedian; + + fprintf (stderr, "norm (1): %f\n", stamps->normValue); + + psFree(stats); + psFree(norms); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate normalization: %f sec", psTimerClear("pmSubtractionCalculateNormalization")); + + return true; +} + +bool pmSubtractionCalculateNormalizationStamp( + pmSubtractionStamp *stamp, // stamp on which to save normalization) + const psKernel *image1, // Input image (target) + const psKernel *image2, // Reference image (convolution source) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ) +{ + double normI1 = 0.0; // Sum of I_1 within the normalisation window (aperture) + double normI2 = 0.0; // Sum of I_2 within the normalisation window (aperture) + double normSquare1 = 0.0; // Sum of (I_1)^2 within the normalisation window (aperture) + double normSquare2 = 0.0; // Sum of (I_2)^2 within the normalisation window (aperture) + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double im1 = image1->kernel[y][x]; + double im2 = image2->kernel[y][x]; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { + normI1 += im1; + normSquare1 += PS_SQR(im1); + } + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { + normI2 += im2; + normSquare2 += PS_SQR(im2); + } + } + } + stamp->norm = normI2 / normI1; + stamp->normSquare1 = normSquare1; + stamp->normSquare2 = normSquare2; + + fprintf (stderr, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); + + return true; +} bool pmSubtractionCalculateEquationThread(psThreadJob *job) @@ -715,5 +592,4 @@ int numKernels = kernels->num; // Number of kernel basis functions int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). @@ -722,5 +598,6 @@ // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the // number of coefficients for the spatial polynomial, normalisation and a constant background offset. - int numParams = numKernels * numSpatial + 1 + numBackground; + // XXX we no longer solve for the normalization and background in the matrix inversion + int numParams = numKernels * numSpatial; if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { // An additional image is convolved @@ -790,18 +667,15 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: - status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, - weight, window, stamp->convolutions1, kernels, - polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image2, stamp->image1, + weight, window, stamp->convolutions1, kernels, polyValues, footprint); break; case PM_SUBTRACTION_MODE_2: - status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, - weight, window, stamp->convolutions2, kernels, - polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode); + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2, + weight, window, stamp->convolutions2, kernels, polyValues, footprint); break; case PM_SUBTRACTION_MODE_DUAL: - status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, - stamp->image1, stamp->image2, + status = calculateDualMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2, weight, window, stamp->convolutions1, stamp->convolutions2, - kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + kernels, polyValues, footprint); break; default: @@ -928,7 +802,8 @@ int numKernels = kernels->num; // Number of kernel basis functions int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for - int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution + // XXX int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + // XXX int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for + int numParams = numKernels * numSpatial; // Number of parameters being solved for + int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { // An additional image is convolved @@ -960,11 +835,10 @@ if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { - // Accumulate the least-squares matricies and vectors + // Accumulate the least-squares matrices and vectors. These are generated for the + // kernel elements, excluding the background and normalization. psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector psVectorInit(sumVector, 0.0); psImageInit(sumMatrix, 0.0); - - psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations int numStamps = 0; // Number of good stamps @@ -976,6 +850,4 @@ (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - psVectorAppend(norms, stamp->norm); - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; @@ -985,108 +857,40 @@ } -#if 0 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); -#endif - psVector *solution = NULL; // Solution to equation! solution = psVectorAlloc(numParams, PS_TYPE_F64); psVectorInit(solution, 0); -#if 0 - // Regular, straight-forward solution - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); -#else - { - // Solve normalisation and background separately - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm - if (!psVectorStats(stats, norms, NULL, NULL, 0)) { - psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); - psFree(stats); - psFree(sumMatrix); - psFree(sumVector); - psFree(norms); - return false; - } - - // double normValue = 1.0; - double normValue = stats->robustMedian; - // double bgValue = 0.0; - - psFree(stats); - -#ifdef TESTING - fprintf(stderr, "Norm: %lf\n", normValue); -#endif - // Solve kernel components - for (int i = 0; i < numSolution1; i++) { - sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; - - sumMatrix->data.F64[i][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i] = 0.0; - } - sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; - sumMatrix->data.F64[bgIndex][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][bgIndex] = 0.0; - - sumMatrix->data.F64[normIndex][normIndex] = 1.0; - sumVector->data.F64[normIndex] = 0.0; - - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - - solution->data.F64[normIndex] = normValue; - } -# endif - -#if (1) + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + for (int i = 0; i < solution->n; i++) { - fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]); - } -#endif + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]); + } if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64); + kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg psVectorInit(kernels->solution1, 0.0); } - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - int spatialOrder = kernels->spatialOrder; // Order of spatial variation - int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms - for (int i = 0; i < numKernels * numPoly; i++) { - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } - - psFree(norms); + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + + // Apply the normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + kernels->solution1->data.F64[normIndex] = stamps->normValue; + kernels->solution1->data.F64[bgIndex] = 0.0; + psFree(solution); psFree(sumVector); psFree(sumMatrix); -#ifdef TESTING - // XXX double-check for NAN in data: - for (int ix = 0; ix < kernels->solution1->n; ix++) { - if (!isfinite(kernels->solution1->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - } else { // Dual convolution solution - - // Accumulation of stamp matrices/vectors + // Accumulate the least-squares matrices and vectors. These are generated for the + // kernel elements, excluding the background and normalization. psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); @@ -1094,7 +898,7 @@ psVectorInit(sumVector, 0.0); - psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations - - int numStamps = 0; // Number of good stamps + int numStamps = 0; // Number of good stamps + double normSquare1 = 0.0; // Sum of (I_1)^2 over stamps + double normSquare2 = 0.0; // Sum of (I_2)^2 over stamps for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -1103,8 +907,11 @@ (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - psVectorAppend(norms, stamp->norm); + normSquare1 += stamp->normSquare1; + normSquare2 += stamp->normSquare2; pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; + } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); } } @@ -1117,11 +924,5 @@ #endif -#if 1 - // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); - - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0); -#endif + calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2); psVector *solution = NULL; // Solution to equation! @@ -1129,87 +930,5 @@ psVectorInit(solution, 0); -#if 0 - // Regular, straight-forward solution - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); -#else - { - // Solve normalisation and background separately - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - -#if 0 - psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); - psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); - - normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; - normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; - normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; - - normVector->data.F64[0] = sumVector->data.F64[normIndex]; - normVector->data.F64[1] = sumVector->data.F64[bgIndex]; - - psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); - - double normValue = normSolution->data.F64[0]; - double bgValue = normSolution->data.F64[1]; - - psFree(normMatrix); - psFree(normVector); - psFree(normSolution); -#endif - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm - if (!psVectorStats(stats, norms, NULL, NULL, 0)) { - psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); - psFree(stats); - psFree(sumMatrix); - psFree(sumVector); - psFree(norms); - return false; - } - - double normValue = stats->robustMedian; - - psFree(stats); - -#ifdef TESTING - fprintf(stderr, "Norm: %lf\n", normValue); -#endif - - // Solve kernel components - for (int i = 0; i < numSolution2; i++) { - sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; - sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; - - sumMatrix->data.F64[i][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i] = 0.0; - - sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; - } - sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; - sumMatrix->data.F64[bgIndex][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][bgIndex] = 0.0; - - sumMatrix->data.F64[normIndex][normIndex] = 1.0; - - sumVector->data.F64[normIndex] = 0.0; - -// save the matrix and vector after the NULLs have been set -#if 0 - psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); - psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); - psVectorWriteFile("sumVector.dat", sumVector); - psFree (save); -#endif - - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6); - // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 3e-4); - // psVectorCopy (solution, sumVector, PS_TYPE_F64); - // psMatrixGJSolve(sumMatrix, solution); - solution->data.F64[normIndex] = normValue; - } -#endif - + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6); #if (1) @@ -1219,11 +938,6 @@ #endif - psFree(sumMatrix); - psFree(sumVector); - - psFree(norms); - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64); + kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); psVectorInit (kernels->solution1, 0.0); } @@ -1233,30 +947,20 @@ } - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - for (int i = 0; i < numKernels * numSpatial; i++) { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; - } - } - - - memcpy(kernels->solution1->data.F64, solution->data.F64, - numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], - numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - + int numKernels = kernels->num; + for (int i = 0; i < numKernels * numSpatial; i++) { + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; + } + + // Apply the normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + kernels->solution1->data.F64[normIndex] = stamps->normValue; + kernels->solution1->data.F64[bgIndex] = 0.0; + + psFree(sumMatrix); + psFree(sumVector); psFree(solution); - } @@ -1927,312 +1631,2 @@ return true; } - - -# if 0 - -#ifdef TESTING - psFitsWriteImageSimple("A.fits", sumMatrix, NULL); - psVectorWriteFile ("B.dat", sumVector); -#endif - -# define SVD_ANALYSIS 0 -# define COEFF_SIG 0.0 -# define SVD_TOL 0.0 - - // Use SVD to determine the kernel coeffs (and validate) - if (SVD_ANALYSIS) { - - // We have sumVector and sumMatrix. we are trying to solve the following equation: - // sumMatrix * x = sumVector. - - // we can use any standard matrix inversion to solve this. However, the basis - // functions in general have substantial correlation, so that the solution may be - // somewhat poorly determined or unstable. If not numerically ill-conditioned, the - // system of equations may be statistically ill-conditioned. Noise in the image - // will drive insignificant, but correlated, terms in the solution. To avoid these - // problems, we can use SVD to identify numerically unconstrained values and to - // avoid statistically badly determined value. - - // A = sumMatrix, B = sumVector - // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T - // x = V (1/w) (U^T B) - // \sigma_x = sqrt(diag(A^{-1})) - // solve for x and A^{-1} to get x & dx - // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 - // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 - - // If I use the SVD trick to re-condition the matrix, I need to break out the - // kernel and normalization terms from the background term. - // XXX is this true? or was this due to an error in the analysis? - - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - - // now pull out the kernel elements into their own square matrix - psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); - psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { - if (ix == bgIndex) continue; - for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { - if (iy == bgIndex) continue; - kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; - ky++; - } - kernelVector->data.F64[kx] = sumVector->data.F64[ix]; - kx++; - } - - psImage *U = NULL; - psImage *V = NULL; - psVector *w = NULL; - if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { - psError(psErrorCodeLast(), false, "failed to perform SVD on sumMatrix\n"); - return NULL; - } - - // calculate A_inverse: - // Ainv = V * w * U^T - psImage *wUt = p_pmSubSolve_wUt (w, U); - psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); - psImage *Xvar = NULL; - psFree (wUt); - -# ifdef TESTING - // kernel terms: - for (int i = 0; i < w->n; i++) { - fprintf (stderr, "w: %f\n", w->data.F64[i]); - } -# endif - // loop over w adding in more and more of the values until chisquare is no longer - // dropping significantly. - // XXX this does not seem to work very well: we seem to need all terms even for - // simple cases... - - psVector *Xsvd = NULL; - { - psVector *Ax = NULL; - psVector *UtB = NULL; - psVector *wUtB = NULL; - - psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); - psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); - psVectorInit (wMask, 1); // start by masking everything - - double chiSquareLast = NAN; - int maxWeight = 0; - - double Axx, Bx, y2; - - // XXX this is an attempt to exclude insignificant modes. - // it was not successful with the ISIS kernel set: removing even - // the least significant mode leaves additional ringing / noise - // because the terms are so coupled. - for (int k = 0; false && (k < w->n); k++) { - - // unmask the k-th weight - wMask->data.U8[k] = 0; - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; - chiSquareLast = chiSquare; - - // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); - if (k && !maxWeight && (deltaChi < 1.0)) { - maxWeight = k; - } - } - - // keep all terms or we get extra ringing - maxWeight = w->n; - psVectorInit (wMask, 1); - for (int k = 0; k < maxWeight; k++) { - wMask->data.U8[k] = 0; - } - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); - - // re-calculate A^{-1} to get new variances: - // Ainv = V * w * U^T - // XXX since we keep all terms, this is identical to Ainv - psImage *wUt = p_pmSubSolve_wUt (wApply, U); - Xvar = p_pmSubSolve_VwUt (V, wUt); - psFree (wUt); - - psFree (Ax); - psFree (UtB); - psFree (wUtB); - psFree (wApply); - psFree (wMask); - } - - // copy the kernel solutions to the full solution vector: - solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { - if (ix == bgIndex) { - solution->data.F64[ix] = 0; - solutionErr->data.F64[ix] = 0.001; - continue; - } - solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); - solution->data.F64[ix] = Xsvd->data.F64[kx]; - kx++; - } - - psFree (kernelMatrix); - psFree (kernelVector); - - psFree (U); - psFree (V); - psFree (w); - - psFree (Ainv); - psFree (Xsvd); - } else { - psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); - if (!luMatrix) { - psError(PM_ERR_DATA, true, "LU Decomposition of least-squares matrix failed.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - psFree(luMatrix); - psFree(permutation); - return NULL; - } - - solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); - psFree(luMatrix); - psFree(permutation); - if (!solution) { - psError(PM_ERR_DATA, true, "Failed to solve the least-squares system.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - return NULL; - } - - // XXX LUD does not provide A^{-1}? fake the error for now - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - for (int ix = 0; ix < sumVector->n; ix++) { - solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; - } - } - - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); - psVectorInit (kernels->solution1, 0.0); - } - - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - int spatialOrder = kernels->spatialOrder; // Order of spatial variation - int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms - for (int i = 0; i < numKernels * numPoly; i++) { - // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); - if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { - // XXX fprintf (stderr, "drop\n"); - kernels->solution1->data.F64[i] = 0.0; - } else { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } - } - // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); - // fprintf (stderr, "chi square Brute = %f\n", chiSquare); - - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); -# endif - -#ifdef TESTING - // XXX double-check for NAN in data: - for (int iy = 0; iy < stamp->matrix->numRows; iy++) { - for (int ix = 0; ix < stamp->matrix->numCols; ix++) { - if (!isfinite(stamp->matrix->data.F64[iy][ix])) { - fprintf (stderr, "WARNING: NAN in matrix\n"); - } - } - } - for (int ix = 0; ix < stamp->vector->n; ix++) { - if (!isfinite(stamp->vector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } - { - psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); - psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); - psFree(inverse); - } - { - psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); - psImage *Xt = psMatrixTranspose(NULL, X); - psImage *XtX = psMatrixMultiply(NULL, Xt, X); - psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); - psFree(X); - psFree(Xt); - psFree(XtX); - } -#endif - Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.h (revision 29164) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.h (revision 29165) @@ -65,3 +65,16 @@ ); +bool pmSubtractionCalculateNormalization( + pmSubtractionStampList *stamps, + const pmSubtractionMode mode); + +bool pmSubtractionCalculateNormalizationStamp( + pmSubtractionStamp *stamp, // stamp on which to save normalization) + const psKernel *input, // Input image (target) + const psKernel *reference, // Reference image (convolution source) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ); + #endif Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c (revision 29164) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c (revision 29165) @@ -28,10 +28,5 @@ static bool useFFT = true; // Do convolutions using FFT -# define SEPARATE 0 -# if (SEPARATE) -# define SUBMODE PM_SUBTRACTION_EQUATION_NORM -# else # define SUBMODE PM_SUBTRACTION_EQUATION_ALL -# endif //#define TESTING @@ -672,5 +667,4 @@ kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, inner, binning, ringsOrder, penalty, bounds, subMode); - // pmSubtractionVisualShowKernels(kernels); } @@ -678,27 +672,4 @@ if (subMode == PM_SUBTRACTION_MODE_UNSURE) { -#if 0 - // Get backgrounds - psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background - psVector *buffer = NULL;// Buffer for stats - if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) { - psError(PM_ERR_DATA, false, "Unable to measure background of image 1."); - psFree(bgStats); - psFree(buffer); - goto MATCH_ERROR; - } - float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1 - if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) { - psError(PM_ERR_DATA, false, "Unable to measure background of image 2."); - psFree(bgStats); - psFree(buffer); - goto MATCH_ERROR; - } - float bg2 = psStatsGetValue(bgStats, BG_STAT); // Background for image 2 - psFree(bgStats); - psFree(buffer); - - pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use -#endif pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); switch (newMode) { @@ -732,6 +703,13 @@ } - // XXX step 1: calculate normalization - psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); + // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue + psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); + if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + // step 1: generate the elements of the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); @@ -739,5 +717,6 @@ } - psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n"); + // step 2: solve the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Solving kernel equations...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); @@ -746,23 +725,4 @@ memCheck(" solve equation"); -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters - psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" calculate equation"); - - psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" solve equation"); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -770,5 +730,4 @@ goto MATCH_ERROR; } - memCheck(" calculate deviations"); @@ -787,5 +746,6 @@ // if we hit the max number of iterations and we have rejected stamps, re-solve if (numRejected > 0) { - // XXX step 1: calculate normalization + + // step 1: generate the elements of the matrix equation Ax = B psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { @@ -794,5 +754,5 @@ } - // solve normalization + // step 2: solve the matrix equation Ax = B psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { @@ -801,23 +761,4 @@ } -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters - psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - // solve kernel parameters - psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" solve equation"); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -837,5 +778,4 @@ goto MATCH_ERROR; } - memCheck("diag outputs"); @@ -1134,21 +1074,4 @@ } -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } -# endif - psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations @@ -1181,22 +1104,4 @@ } psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); - -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations @@ -1301,5 +1206,5 @@ PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); -// float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference + // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference float scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c (revision 29164) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c (revision 29165) @@ -244,4 +244,5 @@ list->flux = NULL; list->normFrac = normFrac; + list->normValue = NAN; list->window1 = NULL; list->window2 = NULL; @@ -343,4 +344,6 @@ stamp->vector = NULL; stamp->norm = NAN; + stamp->normSquare1 = NAN; + stamp->normSquare2 = NAN; return stamp; Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h (revision 29164) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h (revision 29165) @@ -25,8 +25,9 @@ int footprint; ///< Half-size of stamps float normFrac; ///< Fraction of flux in window for normalisation window + float normValue; ///< calculated normalization psKernel *window1; ///< window function generated from ensemble of stamps (input 1) psKernel *window2; ///< window function generated from ensemble of stamps (input 2) - float normWindow1; ///< Size of window for measuring normalisation - float normWindow2; ///< Size of window for measuring normalisation + float normWindow1; ///< Size of window for measuring normalisation + float normWindow2; ///< Size of window for measuring normalisation float sysErr; ///< Systematic error float skyErr; ///< increase effective readnoise @@ -85,4 +86,6 @@ psVector *vector; ///< Least-squares vector, or NULL double norm; ///< Normalisation difference + double normSquare1; ///< Sum(flux^2) for image 1 (used for penalty) + double normSquare2; ///< Sum(flux^2) for image 2 (used for penalty) pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp;