Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26329) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26330) @@ -35,5 +35,6 @@ const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint // (Half-)Size of stamp + int footprint, // (Half-)Size of stamp + const bool normOnly ) { @@ -53,4 +54,5 @@ // Evaluate polynomial-polynomial terms + // XXX we can skip this for normOnly for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { @@ -67,76 +69,89 @@ - for (int i = 0; i < numKernels; i++) { - psKernel *iConv = convolutions->data[i]; // Convolution for index i - for (int j = i; j < numKernels; j++) { - psKernel *jConv = convolutions->data[j]; // Convolution for index j - - double sumCC = 0.0; // Sum of convolution products - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; + // if we only want to calculate the normalization (and background model?) + // then we should assign these values in the matrix and vector to 1 (diagonal) or 0 (off-diagonal) + // XXX need to disable normalization calculation if !normOnly + if (normOnly) { + psImageInit(matrix, 0.0); + psVectorInit(vector, 1.0); + for (int i = 0; i < matrix->numCols; i++) { + matrix->data.F64[i][i] = 1.0; + } + } else { + for (int i = 0; i < numKernels; i++) { + psKernel *iConv = convolutions->data[i]; // Convolution for index i + for (int j = i; j < numKernels; j++) { + psKernel *jConv = convolutions->data[j]; // Convolution for index j + + double sumCC = 0.0; // Sum of convolution products + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; #ifdef USE_WEIGHT - cc *= weight->kernel[y][x]; + cc *= weight->kernel[y][x]; +#endif +#ifdef USE_WINDOW + if (window) { + cc *= window->kernel[y][x]; + } +#endif + sumCC += cc; + } + } + + // Spatial variation + 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++) { + float conv = iConv->kernel[y][x]; + float in = input->kernel[y][x]; + float ref = reference->kernel[y][x]; + double ic = in * conv; + double rc = ref * conv; + double c = conv; +#ifdef USE_WEIGHT + float wtVal = weight->kernel[y][x]; + ic *= wtVal; + rc *= wtVal; + c *= wtVal; #endif #ifdef USE_WINDOW if (window) { - cc *= window->kernel[y][x]; + float winVal = window->kernel[y][x]; + ic *= winVal; + rc *= winVal; + c *= winVal; } #endif - sumCC += cc; - } - } - - // Spatial variation - 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++) { - float conv = iConv->kernel[y][x]; - float in = input->kernel[y][x]; - float ref = reference->kernel[y][x]; - double ic = in * conv; - double rc = ref * conv; - double c = conv; -#ifdef USE_WEIGHT - float wtVal = weight->kernel[y][x]; - ic *= wtVal; - rc *= wtVal; - c *= wtVal; -#endif -#ifdef USE_WINDOW - if (window) { - float winVal = window->kernel[y][x]; - ic *= winVal; - rc *= winVal; - c *= winVal; + sumIC += ic; + sumRC += rc; + sumC += c; } -#endif - 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]; - matrix->data.F64[iIndex][normIndex] = normTerm; - matrix->data.F64[normIndex][iIndex] = normTerm; - matrix->data.F64[iIndex][bgIndex] = bgTerm; - matrix->data.F64[bgIndex][iIndex] = bgTerm; - vector->data.F64[iIndex] = sumIC * poly[iTerm]; - } - } + } + // Spatial variation + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + double normTerm = sumRC * poly[iTerm]; + double bgTerm = sumC * poly[iTerm]; + matrix->data.F64[iIndex][normIndex] = normTerm; + matrix->data.F64[normIndex][iIndex] = normTerm; + matrix->data.F64[iIndex][bgIndex] = bgTerm; + matrix->data.F64[bgIndex][iIndex] = bgTerm; + vector->data.F64[iIndex] = sumIC * poly[iTerm]; + } + } + } + + // XXX need to disable normalization calculation if !normOnly double sumRR = 0.0; // Sum of the reference product @@ -574,9 +589,9 @@ int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index - return pmSubtractionCalculateEquationStamp(stamps, kernels, index); + return pmSubtractionCalculateEquationStamp(stamps, kernels, index, normOnly); } bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, - int index) + int index, bool normOnly) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -591,4 +606,7 @@ 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). + // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3 + // 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. @@ -598,5 +616,5 @@ psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state."); - // Generate convolutions + // Generate convolutions: these are generated once and saved if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", index); @@ -648,12 +666,13 @@ status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, stamp->weight, stamps->window, stamp->convolutions1, kernels, - polyValues, footprint); + polyValues, footprint, normOnly); break; case PM_SUBTRACTION_MODE_2: status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2, stamp->weight, stamps->window, stamp->convolutions2, kernels, - polyValues, footprint); + polyValues, footprint, normOnly); break; case PM_SUBTRACTION_MODE_DUAL: + psAbort ("dual is disabled: need to add window and normOnly"); if (new) { stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); @@ -739,5 +758,5 @@ } -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels) +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const bool normOnly) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -759,4 +778,5 @@ psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, normOnly, PS_TYPE_BOOL); if (!psThreadJobAddPending(job)) { psFree(job); @@ -765,5 +785,5 @@ psFree(job); } else { - pmSubtractionCalculateEquationStamp(stamps, kernels, i); + pmSubtractionCalculateEquationStamp(stamps, kernels, i, normOnly); } } Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h (revision 26329) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h (revision 26330) @@ -12,10 +12,12 @@ bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps const pmSubtractionKernels *kernels, ///< Kernel parameters - int index ///< Index of stamp + int index, ///< Index of stamp + const bool normOnly ); /// Calculate the least-squares equation to match the image quality bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps - const pmSubtractionKernels *kernels ///< Kernel parameters + const pmSubtractionKernels *kernels, ///< Kernel parameters + const bool normOnly ); Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26329) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26330) @@ -560,4 +560,13 @@ } + // XXX step 1: calculate normalization + psTrace("psModules.imcombine", 3, "Calculating equation...\n"); + if (!pmSubtractionCalculateEquation(stamps, kernels, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + // XXX step 2: calculate only other parameters +# if (0) psTrace("psModules.imcombine", 3, "Calculating equation...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels)) { @@ -565,4 +574,5 @@ goto MATCH_ERROR; } +# endif memCheck(" calculate equation"); Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionThreads.c (revision 26329) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionThreads.c (revision 26330) @@ -69,5 +69,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 3); + psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 4); task->function = &pmSubtractionCalculateEquationThread; psThreadTaskAdd(task);