Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26331) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26332) @@ -462,5 +462,5 @@ { psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data -#if 1 +#if 0 // Convolving using separable convolution psVector *xKernel = preCalc->data[0]; // Kernel in x Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26331) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26332) @@ -26,4 +26,5 @@ // Calculate the least-squares matrix and vector +// XXX why are we calculating these values on each iteration? aren't they identical (no change in pixel values...) static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated psVector *vector, // Least-squares vector, updated @@ -36,5 +37,5 @@ const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp - const bool normOnly + const pmSubtractionEquationCalculationMode mode ) { @@ -54,5 +55,5 @@ // Evaluate polynomial-polynomial terms - // XXX we can skip this for normOnly + // 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++) { @@ -68,37 +69,40 @@ } - - // 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 + // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we + // choose to calculate + psImageInit(matrix, 0.0); + psVectorInit(vector, 1.0); + for (int i = 0; i < matrix->numCols; i++) { + 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] + // normalization + // bg 0, bg 1, bg 2 (only 0 is currently used?) + + 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 (weight) { cc *= weight->kernel[y][x]; -#endif -#ifdef USE_WINDOW - if (window) { - cc *= window->kernel[y][x]; - } -#endif - sumCC += cc; } + if (window) { + cc *= window->kernel[y][x]; + } + sumCC += cc; } - - // Spatial variation + } + + // 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) { @@ -109,49 +113,51 @@ } } - - 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 + } + + 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; + if (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; - } -#endif - sumIC += ic; - sumRC += rc; - sumC += c; } + if (window) { + float winVal = window->kernel[y][x]; + 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]; + } + // 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]; } } } - - // XXX need to disable normalization calculation if !normOnly double sumRR = 0.0; // Sum of the reference product @@ -167,13 +173,12 @@ double rr = PS_SQR(ref); double one = 1.0; -#ifdef USE_WEIGHT - float wtVal = weight->kernel[y][x]; - rr *= wtVal; - ir *= wtVal; - in *= wtVal; - ref *= wtVal; - one *= wtVal; -#endif -#ifdef USE_WINDOW + 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]; @@ -184,5 +189,4 @@ one *= winVal; } -#endif sumRR += rr; sumIR += ir; @@ -192,10 +196,16 @@ } } - matrix->data.F64[normIndex][normIndex] = sumRR; - matrix->data.F64[bgIndex][bgIndex] = sum1; - matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR; - vector->data.F64[normIndex] = sumIR; - vector->data.F64[bgIndex] = sumI; - + if (mode | PM_SUBTRACTION_EQUATION_NORM) { + matrix->data.F64[normIndex][normIndex] = sumRR; + vector->data.F64[normIndex] = sumIR; + } + 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; + } return true; } @@ -267,10 +277,10 @@ double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; -#ifdef USE_WEIGHT - float wtVal = weight->kernel[y][x]; - aa *= wtVal; - bb *= wtVal; - ab *= wtVal; -#endif + if (weight) { + float wtVal = weight->kernel[y][x]; + aa *= wtVal; + bb *= wtVal; + ab *= wtVal; + } sumAA += aa; sumBB += bb; @@ -299,7 +309,7 @@ for (int x = - footprint; x <= footprint; x++) { double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; -#ifdef USE_WEIGHT - ab *= weight->kernel[y][x]; -#endif + if (weight) { + ab *= weight->kernel[y][x]; + } sumAB += ab; } @@ -336,16 +346,15 @@ double i1i2 = i1 * i2; -#ifdef USE_WEIGHT - float wtVal = weight->kernel[y][x]; - ai2 *= wtVal; - bi2 *= wtVal; - ai1 *= wtVal; - bi1 *= wtVal; - i1i2 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; -#endif - + if (weight) { + float wtVal = weight->kernel[y][x]; + ai2 *= wtVal; + bi2 *= wtVal; + ai1 *= wtVal; + bi1 *= wtVal; + i1i2 *= wtVal; + a *= wtVal; + b *= wtVal; + i2 *= wtVal; + } sumAI2 += ai2; sumBI2 += bi2; @@ -392,13 +401,12 @@ double i1i2 = i1 * i2; -#ifdef USE_WEIGHT - float wtVal = weight->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; -#endif - + if (weight) { + float wtVal = weight->kernel[y][x]; + i1 *= wtVal; + i1i1 *= wtVal; + one *= wtVal; + i2 *= wtVal; + i1i2 *= wtVal; + } sumI1 += i1; sumI1I1 += i1i1; @@ -588,10 +596,11 @@ const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index - - return pmSubtractionCalculateEquationStamp(stamps, kernels, index, normOnly); + pmSubtractionEquationCalculationMode mode = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model + + return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode); } bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, - int index, bool normOnly) + int index, const pmSubtractionEquationCalculationMode mode) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -662,17 +671,28 @@ bool status; // Status of least-squares matrix/vector calculation + + psKernel *weight = NULL; + psKernel *window = NULL; + +#ifdef USE_WEIGHT + weight = stamp->weight; +#endif +#ifdef USE_WINDOW + window = stamps->window; +#endif + switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, - stamp->weight, stamps->window, stamp->convolutions1, kernels, - polyValues, footprint, normOnly); + weight, window, stamp->convolutions1, kernels, + polyValues, footprint, mode); 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, normOnly); + weight, window, stamp->convolutions2, kernels, + polyValues, footprint, mode); break; case PM_SUBTRACTION_MODE_DUAL: - psAbort ("dual is disabled: need to add window and normOnly"); + psAbort ("dual is disabled: need to add window and calculation mode"); if (new) { stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); @@ -686,5 +706,5 @@ #endif status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, - stamp->matrixX, stamp->image1, stamp->image2, stamp->weight, + stamp->matrixX, stamp->image1, stamp->image2, weight, stamp->convolutions1, stamp->convolutions2, kernels, polyValues, footprint); @@ -758,5 +778,5 @@ } -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const bool normOnly) +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -778,5 +798,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); + PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32); if (!psThreadJobAddPending(job)) { psFree(job); @@ -785,5 +805,5 @@ psFree(job); } else { - pmSubtractionCalculateEquationStamp(stamps, kernels, i, normOnly); + pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode); } } @@ -803,5 +823,7 @@ } -bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, + const pmSubtractionStampList *stamps, + const pmSubtractionEquationCalculationMode mode) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); @@ -941,5 +963,37 @@ return NULL; } - kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation); + + psVector *solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + if (!solution) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + return NULL; + } + + 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++) { + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + psFree (solution); #ifdef TESTING @@ -952,11 +1006,4 @@ #endif - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - if (!kernels->solution1) { - psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); - return NULL; - } } else { // Dual convolution solution Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h (revision 26331) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h (revision 26332) @@ -4,4 +4,12 @@ #include "pmSubtractionStamps.h" #include "pmSubtractionKernels.h" + +typedef enum { + PM_SUBTRACTION_EQUATION_NONE = 0x00, + PM_SUBTRACTION_EQUATION_NORM = 0x01, + PM_SUBTRACTION_EQUATION_BG = 0x02, + PM_SUBTRACTION_EQUATION_KERNELS = 0x04, + PM_SUBTRACTION_EQUATION_ALL = 0x05, // value should be NORM | BG | KERNELS +} pmSubtractionEquationCalculationMode; /// Execute a thread job to calculate the least-squares equation for a stamp @@ -13,16 +21,17 @@ const pmSubtractionKernels *kernels, ///< Kernel parameters int index, ///< Index of stamp - const bool normOnly - ); + const pmSubtractionEquationCalculationMode mode + ); /// Calculate the least-squares equation to match the image quality bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps const pmSubtractionKernels *kernels, ///< Kernel parameters - const bool normOnly - ); + const pmSubtractionEquationCalculationMode mode + ); /// Solve the least-squares equation to match the image quality bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters - const pmSubtractionStampList *stamps ///< Stamps + const pmSubtractionStampList *stamps, ///< Stamps + const pmSubtractionEquationCalculationMode mode ); Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26331) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26332) @@ -561,29 +561,35 @@ // XXX step 1: calculate normalization - psTrace("psModules.imcombine", 3, "Calculating equation...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, true)) { + psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); + // if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_NORM)) { + if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); goto MATCH_ERROR; } - // XXX step 2: calculate only other parameters + psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n"); + // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_NORM)) { + if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + memCheck(" solve equation"); + # if (0) - psTrace("psModules.imcombine", 3, "Calculating equation...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels)) { + // XXX step 2: calculate kernel parameters + psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n"); + if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { psError(PS_ERR_UNKNOWN, 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(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + memCheck(" solve equation"); # endif - - memCheck(" calculate equation"); - - psTrace("psModules.imcombine", 3, "Solving equation...\n"); - - if (!pmSubtractionSolveEquation(kernels, stamps)) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - memCheck(" solve equation"); psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations @@ -609,9 +615,13 @@ // if we hit the max number of iterations and we have rejected stamps, re-solve if (numRejected > 0) { - psTrace("psModules.imcombine", 3, "Solving equation...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps)) { + // calculate kernel parameters + psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); + // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { + if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); goto MATCH_ERROR; } + memCheck(" solve equation"); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -918,15 +928,31 @@ assert(kernels); - psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); + // if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_NORM)) { + if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); return false; } - psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps)) { + psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description); + // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_NORM)) { + if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); return false; } + +# if (0) + psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description); + if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { + psError(PS_ERR_UNKNOWN, 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(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } +# endif psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); @@ -949,9 +975,11 @@ // Allow re-fit with reduced stamps set psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps)) { + // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { + if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); return false; } psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26331) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26332) @@ -263,4 +263,5 @@ static psImage *targetImage = NULL; static psImage *residualImage = NULL; +static psImage *fresidualImage = NULL; static psImage *differenceImage = NULL; static psImage *convolutionImage = NULL; @@ -289,4 +290,5 @@ targetImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); residualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + fresidualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); differenceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); @@ -295,4 +297,5 @@ psImageInit (targetImage, 0.0); psImageInit (residualImage, 0.0); + psImageInit (fresidualImage, 0.0); psImageInit (differenceImage, 0.0); psImageInit (convolutionImage, 0.0); @@ -343,4 +346,11 @@ for (int x = -footprint; x <= footprint; x++) { residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm + convolution->kernel[y][x]; + } + } + + // insert the (fresidual) kernel into the (fresidual) image: + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / target->kernel[y][x]; } } @@ -358,6 +368,7 @@ pmVisualScaleImage(kapa1, convolutionImage, "Convolution Stamps", 2, true); + // pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true); + pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true); pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true); - pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true); pmVisualAskUser(NULL); @@ -367,4 +378,5 @@ psFree(differenceImage); psFree(residualImage); + psFree(fresidualImage); targetImage = NULL; @@ -373,4 +385,5 @@ differenceImage = NULL; residualImage = NULL; + fresidualImage = NULL; return true;