Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.c (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.c (revision 28866) @@ -33,4 +33,9 @@ #define USE_KERNEL_ERR // Use kernel error image? #define NUM_COVAR_POS 5 // Number of positions for covariance calculation + +// XXX we need to pass these fwhm values elsewhere. These should go on one of the structure, but +// things are too confusing to do that now. just save them here. +static float FWHM1 = NAN; +static float FWHM2 = NAN; ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -752,5 +757,5 @@ -bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint) +bool pmSubtractionConvolveStamp (pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int footprint) { PS_ASSERT_PTR_NON_NULL(stamp, false); @@ -774,4 +779,7 @@ stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint); + if (!pmSubtractionKernelPenaltiesStamp(stamp, kernels)) { + psAbort("failure in penalties"); + } break; default: @@ -1413,2 +1421,16 @@ return true; } + +bool pmSubtractionGetFWHMs(float *fwhm1, float *fwhm2) { + + *fwhm1 = FWHM1; + *fwhm2 = FWHM2; + return true; +} + +bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2) { + + FWHM1 = fwhm1; + FWHM2 = fwhm2; + return true; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.h (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.h (revision 28866) @@ -59,5 +59,5 @@ /// Convolve the reference stamp with the kernel components bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve - const pmSubtractionKernels *kernels, ///< Kernel parameters + pmSubtractionKernels *kernels, ///< Kernel parameters int footprint ///< Half-size of region over which to calculate equation ); @@ -157,4 +157,7 @@ ); +bool pmSubtractionGetFWHMs(float *fwhm1, float *fwhm2); +bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2); + /// @} #endif Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c (revision 28866) @@ -38,5 +38,6 @@ const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp - int normWindow, // Window (half-)size for normalisation measurement + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2, // Window (half-)size for normalisation measurement const pmSubtractionEquationCalculationMode mode ) @@ -184,6 +185,8 @@ double one = 1.0; - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { normI1 += ref; + } + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { normI2 += in; } @@ -214,4 +217,6 @@ *norm = normI2 / normI1; + + fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); if (mode & PM_SUBTRACTION_EQUATION_NORM) { @@ -262,5 +267,6 @@ const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp - int normWindow, // Window (half-)size for normalisation measurement + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2, // Window (half-)size for normalisation measurement const pmSubtractionEquationCalculationMode mode ) @@ -492,6 +498,8 @@ double i1i2 = i1 * i2; - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { normI1 += i1; + } + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { normI2 += i2; } @@ -522,4 +530,5 @@ *norm = normI2 / normI1; + fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); if (mode & PM_SUBTRACTION_EQUATION_NORM) { @@ -559,8 +568,8 @@ // 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 - ) + psVector *vector, // Vector to which to add in penalty term + const pmSubtractionKernels *kernels, // Kernel parameters + float norm // Normalisation + ) { if (kernels->penalty == 0.0) { @@ -568,5 +577,7 @@ } - psVector *penalties = kernels->penalties; // Penalties for each kernel component + psVector *penalties1 = kernels->penalties1; // Penalties for each kernel component (input) + psVector *penalties2 = kernels->penalties2; // Penalties for each kernel component (ref) + int spatialOrder = kernels->spatialOrder; // Order of spatial variations int numKernels = kernels->num; // Number of kernel components @@ -588,8 +599,10 @@ for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { // Contribution to chi^2: a_i^2 P_i - psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty"); - matrix->data.F64[index][index] += norm * penalties->data.F32[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) { - matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i]; + fprintf (stderr, "penalty: %f + %f (%f * %f)\n", 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 @@ -682,5 +695,5 @@ pmSubtractionStampList *stamps = job->args->data[0]; // List of stamps - const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels + pmSubtractionKernels *kernels = job->args->data[1]; // Kernels int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index pmSubtractionEquationCalculationMode mode = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model @@ -689,5 +702,5 @@ } -bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, +bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, int index, const pmSubtractionEquationCalculationMode mode) { @@ -778,10 +791,10 @@ status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, weight, window, stamp->convolutions1, kernels, - polyValues, footprint, stamps->normWindow, mode); + polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); 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->normWindow, mode); + polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode); break; case PM_SUBTRACTION_MODE_DUAL: @@ -789,5 +802,5 @@ stamp->image1, stamp->image2, weight, window, stamp->convolutions1, stamp->convolutions2, - kernels, polyValues, footprint, stamps->normWindow, mode); + kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); break; default: @@ -830,5 +843,5 @@ } -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode) { @@ -996,4 +1009,5 @@ } + // double normValue = 1.0; double normValue = stats->robustMedian; // double bgValue = 0.0; @@ -1023,4 +1037,10 @@ } # endif + +#if (1) + for (int i = 0; i < solution->n; i++) { + fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]); + } +#endif if (!kernels->solution1) { @@ -1096,5 +1116,5 @@ int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0); + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 10000.0); #endif @@ -1177,5 +1197,5 @@ -#ifdef TESTING +#if (1) for (int i = 0; i < solution->n; i++) { fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.h (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.h (revision 28866) @@ -19,5 +19,5 @@ /// Calculate the least-squares equation to match the image quality for a single stamp bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps - const pmSubtractionKernels *kernels, ///< Kernel parameters + pmSubtractionKernels *kernels, ///< Kernel parameters int index, ///< Index of stamp const pmSubtractionEquationCalculationMode mode @@ -26,5 +26,5 @@ /// Calculate the least-squares equation to match the image quality bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps - const pmSubtractionKernels *kernels, ///< Kernel parameters + pmSubtractionKernels *kernels, ///< Kernel parameters const pmSubtractionEquationCalculationMode mode ); Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c (revision 28866) @@ -26,5 +26,6 @@ psFree(kernels->vStop); psFree(kernels->preCalc); - psFree(kernels->penalties); + psFree(kernels->penalties1); + psFree(kernels->penalties2); psFree(kernels->solution1); psFree(kernels->solution2); @@ -140,5 +141,8 @@ kernels->v = psVectorRealloc(kernels->v, start + numNew); kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew); - kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew); + + kernels->penalties1 = psVectorRealloc(kernels->penalties1, start + numNew); + kernels->penalties2 = psVectorRealloc(kernels->penalties2, start + numNew); + kernels->inner = start; kernels->num += numNew; @@ -156,6 +160,12 @@ kernels->v->data.S32[index] = v; kernels->preCalc->data[index] = NULL; - kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); - psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty"); + + // XXX this needs to be changed to use the *convolved* second moment + kernels->penalties1->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); + psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty"); + + kernels->penalties2->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); + psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty"); + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); } @@ -166,12 +176,14 @@ kernels->v->n = start + numNew; kernels->preCalc->n = start + numNew; - kernels->penalties->n = start + numNew; + + kernels->penalties1->n = start + numNew; + kernels->penalties2->n = start + numNew; return true; } -bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, - int index, int size, int uOrder, int vOrder, float fwhm, - bool AlardLuptonStyle, bool forceZeroNull) +static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, + int index, int uOrder, int vOrder, float fwhm, + bool AlardLuptonStyle, bool forceZeroNull) { // we have 4 cases here: @@ -182,14 +194,13 @@ // Calculate moments - double penalty = 0.0; // Moment, for penalty double sum = 0.0, sum2 = 0.0; // Sum of kernel component float min = INFINITY, max = -INFINITY; // Minimum and maximum kernel value - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { + + for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) { + for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) { double value = preCalc->kernel->kernel[v][u]; double value2 = PS_SQR(value); sum += value; sum2 += value2; - penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v))); min = PS_MIN(value, min); max = PS_MAX(value, max); @@ -198,5 +209,5 @@ #if 0 - fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty); + fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max); #endif @@ -239,5 +250,4 @@ psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32)); - penalty *= 1.0 / sum2; if (zeroNull) { @@ -249,6 +259,6 @@ double sum = 0.0; // Sum of kernel component float min = INFINITY, max = -INFINITY; // Minimum and maximum kernel value - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { + for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) { + for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) { sum += preCalc->kernel->kernel[v][u]; min = PS_MIN(preCalc->kernel->kernel[v][u], min); @@ -267,7 +277,5 @@ } kernels->preCalc->data[index] = preCalc; - kernels->penalties->data.F32[index] = kernels->penalty * penalty; - psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty"); - psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty); + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder); return true; @@ -321,6 +329,6 @@ pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); - // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); + // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false); } } @@ -379,5 +387,5 @@ for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); } } @@ -385,5 +393,5 @@ // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian) pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true); } } @@ -437,5 +445,5 @@ for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); } } @@ -506,5 +514,5 @@ // XXX do we use Alard-Lupton normalization (last param true) or not? - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); // XXXX test demo that deconvolved kernel is valid @@ -572,5 +580,9 @@ kernels->preCalc = psArrayAlloc(numBasisFunctions); kernels->penalty = penalty; - kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); + kernels->penalties1 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); + psVectorInit(kernels->penalties1, NAN); + kernels->penalties2 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); + psVectorInit(kernels->penalties2, NAN); + kernels->havePenalties = false; kernels->size = size; kernels->inner = 0; @@ -771,5 +783,6 @@ psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels."); - psVectorInit(kernels->penalties, 0.0); + psVectorInit(kernels->penalties1, 0.0); + psVectorInit(kernels->penalties2, 0.0); return kernels; @@ -866,5 +879,6 @@ psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels."); - psVectorInit(kernels->penalties, 0.0); + psVectorInit(kernels->penalties1, 0.0); + psVectorInit(kernels->penalties2, 0.0); return kernels; @@ -1040,6 +1054,12 @@ kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; - kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); - if (!isfinite(kernels->penalties->data.F32[index])) { + + // XXX convert to use the convolved 2nd moment + kernels->penalties1->data.F32[index] = kernels->penalty * fabsf(moment); + if (!isfinite(kernels->penalties1->data.F32[index])) { + psAbort ("invalid penalty"); + } + kernels->penalties2->data.F32[index] = kernels->penalty * fabsf(moment); + if (!isfinite(kernels->penalties2->data.F32[index])) { psAbort ("invalid penalty"); } @@ -1247,5 +1267,6 @@ out->preCalc = psMemIncrRefCounter(in->preCalc); out->penalty = in->penalty; - out->penalties = psMemIncrRefCounter(in->penalties); + out->penalties1 = psMemIncrRefCounter(in->penalties1); + out->penalties2 = psMemIncrRefCounter(in->penalties2); out->uStop = psMemIncrRefCounter(in->uStop); out->vStop = psMemIncrRefCounter(in->vStop); Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.h (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.h (revision 28866) @@ -39,5 +39,7 @@ psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM) float penalty; ///< Penalty for wideness - psVector *penalties; ///< Penalty for each kernel component + psVector *penalties1; ///< Penalty for each kernel component + psVector *penalties2; ///< Penalty for each kernel component + bool havePenalties; ///< flag to test if we have already calculated the penalties or not. int size; ///< The half-size of the kernel int inner; ///< The size of an inner region @@ -308,4 +310,3 @@ ); - #endif Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.c (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.c (revision 28866) @@ -1304,4 +1304,7 @@ float scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor + // XXX save these values in a static for later use + pmSubtractionSetFWHMs(fwhm1, fwhm2); + if (isfinite(scaleMin) && scale < scaleMin) { scale = scaleMin; Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.h (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.h (revision 28866) @@ -110,4 +110,3 @@ ); - #endif Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c (revision 28866) @@ -51,5 +51,6 @@ psFree(list->y); psFree(list->flux); - psFree(list->window); + psFree(list->window1); + psFree(list->window2); } @@ -230,7 +231,9 @@ list->y = NULL; list->flux = NULL; - list->window = NULL; list->normFrac = normFrac; - list->normWindow = 0; + list->window1 = NULL; + list->window2 = NULL; + list->normWindow1 = 0; + list->normWindow2 = 0; list->footprint = footprint; list->sysErr = sysErr; @@ -253,7 +256,9 @@ out->y = NULL; out->flux = NULL; - out->window = psMemIncrRefCounter(in->window); + out->window1 = psMemIncrRefCounter(in->window1); + out->window2 = psMemIncrRefCounter(in->window2); out->footprint = in->footprint; - out->normWindow = in->normWindow; + out->normWindow1 = in->normWindow1; + out->normWindow2 = in->normWindow2; for (int i = 0; i < num; i++) { @@ -643,7 +648,11 @@ int size = stamps->footprint; // Size of postage stamps - psFree (stamps->window); - stamps->window = psKernelAlloc(-size, size, -size, size); - psImageInit(stamps->window->image, 0.0); + psFree (stamps->window1); + stamps->window1 = psKernelAlloc(-size, size, -size, size); + psImageInit(stamps->window1->image, 0.0); + + psFree (stamps->window2); + stamps->window2 = psKernelAlloc(-size, size, -size, size); + psImageInit(stamps->window2->image, 0.0); // generate normalizations for each stamp @@ -674,6 +683,8 @@ // generate the window pixels - double sum = 0.0; // Sum inside the window - float maxValue = 0.0; // Maximum value, for normalisation + double sum1 = 0.0; // Sum inside the window + double sum2 = 0.0; // Sum inside the window + float maxValue1 = 0.0; // Maximum value, for normalisation + float maxValue2 = 0.0; // Maximum value, for normalisation for (int y = -size; y <= size; y++) { for (int x = -size; x <= size; x++) { @@ -696,4 +707,5 @@ } float f1 = stats->robustMedian; + psStatsInit (stats); if (!psVectorStats (stats, flux2, NULL, NULL, 0)) { @@ -702,34 +714,104 @@ float f2 = stats->robustMedian; - stamps->window->kernel[y][x] = f1 + f2; + stamps->window1->kernel[y][x] = f1; if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) { - sum += stamps->window->kernel[y][x]; - } - maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]); - } - } - - psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n", - sum, (1.0 - stamps->normFrac) * sum); - bool done = false; - for (int radius = 1; radius <= size && !done; radius++) { - double within = 0.0; + sum1 += stamps->window1->kernel[y][x]; + } + maxValue1 = PS_MAX(maxValue1, stamps->window1->kernel[y][x]); + + stamps->window2->kernel[y][x] = f2; + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) { + sum2 += stamps->window2->kernel[y][x]; + } + maxValue2 = PS_MAX(maxValue2, stamps->window2->kernel[y][x]); + } + } + + psTrace("psModules.imcombine", 3, "Window total (1): %f, threshold: %f\n", sum1, (1.0 - stamps->normFrac) * sum1); + psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2); + +# if (0) + // this block attempts to calculate the radius based on the first radial moment + bool done1 = false; + bool done2 = false; + double prior1 = 0.0; + double prior2 = 0.0; + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + float r = hypot(x, y); + Sr1 += r * stamps->window1->kernel[y][x]; + Sr2 += r * stamps->window2->kernel[y][x]; + Sf1 += stamps->window1->kernel[y][x]; + Sf2 += stamps->window2->kernel[y][x]; + } + } + + float R1 = Sr1 / Sf1; + float R2 = Sr2 / Sf2; + + stamps->normWindow1 = 2.5*R1; + stamps->normWindow1 = 2.5*R2; +# else + // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent). + // It did not do very well (though a true curve-of-growth analysis might be better...) + bool done1 = false; + bool done2 = false; + double prior1 = 0.0; + double prior2 = 0.0; + double delta1o = 1.0; + double delta2o = 1.0; + for (int radius = 1; radius <= size && !(done1 && done2); radius++) { + double within1 = 0.0; + double within2 = 0.0; for (int y = -radius; y <= radius; y++) { for (int x = -radius; x <= radius; x++) { if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { - within += stamps->window->kernel[y][x]; + within1 += stamps->window1->kernel[y][x]; } - } - } - psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within); - if (within > (1.0 - stamps->normFrac) * sum) { - stamps->normWindow = radius; - done = true; - } - } - - psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow); - if (stamps->normWindow == 0 || stamps->normWindow >= size) { - psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size."); + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { + within2 += stamps->window2->kernel[y][x]; + } + } + } + double delta1 = (within1 - prior1) / within1; + if (!done1 && (fabs(delta1) < stamps->normFrac)) { + // interpolate to the radius at which delta2 is normFrac: + stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1); + fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1); + done1 = true; + } + double delta2 = (within2 - prior2) / within2; + if (!done2 && (fabs(delta2) < stamps->normFrac)) { + // interpolate to the radius at which delta2 is normFrac: + stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2); + fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2); + done2 = true; + } + psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2); + + prior1 = within1; + prior2 = within2; + delta1o = delta1; + delta2o = delta2; + + // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) { + // stamps->normWindow1 = radius; + // done1 = true; + // } + // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) { + // stamps->normWindow2 = radius; + // done2 = true; + // } + + } +# endif + + psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2); + if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) { + psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1)."); + return false; + } + if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) { + psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2)."); return false; } @@ -738,12 +820,22 @@ for (int y = -size; y <= size; y++) { for (int x = -size; x <= size; x++) { - stamps->window->kernel[y][x] /= maxValue; - } - } - -#if 0 + stamps->window1->kernel[y][x] /= maxValue1; + } + } + // re-normalize so chisquare values are sensible + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + stamps->window2->kernel[y][x] /= maxValue2; + } + } + +#if 1 { - psFits *fits = psFitsOpen ("window.fits", "w"); - psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL); + psFits *fits = NULL; + fits = psFitsOpen ("window1.fits", "w"); + psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL); + psFitsClose (fits); + fits = psFitsOpen ("window2.fits", "w"); + psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL); psFitsClose (fits); } @@ -752,8 +844,100 @@ psFree (stats); psFree (flux1); - psFree(flux2); + psFree (flux2); psFree (norm1); psFree (norm2); return true; +} + +static pthread_mutex_t getPenaltiesMutex = PTHREAD_MUTEX_INITIALIZER; + +// kernels->penalty is an overall scaling factor (user-supplied) +bool pmSubtractionKernelPenaltiesStamp(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels) +{ + // we only need the penalties if we are doing dual convolution + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true; + + // we only calculate the penalties once. + if (kernels->havePenalties) return true; + + // in a threaded context, only one thread can calculate the penalties. attempt to grab a + // mutex before continuing + pthread_mutex_lock(&getPenaltiesMutex); + + // did someone else already get the mutex and do this? + if (kernels->havePenalties) { + pthread_mutex_unlock(&getPenaltiesMutex); + return true; + } + + for (int i = 0; i < kernels->num; i++) { + pmSubtractionKernelPenalties(stamp, kernels, i); + } + + kernels->havePenalties = true; + pthread_mutex_unlock(&getPenaltiesMutex); + return true; +} + +# define EMPIRICAL 0 + +// kernels->penalty is an overall scaling factor (user-supplied) +bool pmSubtractionKernelPenalties(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int index) +{ + float penalty1, penalty2; + float fwhm1, fwhm2; + + // XXX this is annoyingly hack-ish + pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + + if (EMPIRICAL) { + psKernel *convolution1 = stamp->convolutions1->data[index]; + penalty1 = pmSubtractionKernelPenaltySingle(convolution1); + + psKernel *convolution2 = stamp->convolutions2->data[index]; + penalty2 = pmSubtractionKernelPenaltySingle(convolution2); + } else { + pmSubtractionKernelPreCalc *kernel = kernels->preCalc->data[index]; + float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel); + + penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment + penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment + } + kernels->penalties1->data.F32[index] = kernels->penalty * penalty1; + psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty"); + + kernels->penalties2->data.F32[index] = kernels->penalty * penalty2; + psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty"); + + fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2); + + return true; +} + +float pmSubtractionKernelPenaltySingle(psKernel *kernel) +{ + // Calculate moments + double penalty = 0.0; // Moment, for penalty + double sum = 0.0, sum2 = 0.0; // Sum of kernel component + float min = INFINITY, max = -INFINITY; // Minimum and maximum kernel value + for (int v = kernel->yMin; v <= kernel->yMax; v++) { + for (int u = kernel->xMin; u <= kernel->xMax; u++) { + double value = kernel->kernel[v][u]; + double value2 = PS_SQR(value); + sum += value; + sum2 += value2; + penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v))); + min = PS_MIN(value, min); + max = PS_MAX(value, max); + } + } + penalty *= 1.0 / sum2; + + if (0) { + fprintf(stderr, "min: %lf, max: %lf, moment: %lf\n", min, max, penalty); + // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty); + } + + return penalty; } Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h (revision 28866) @@ -24,7 +24,9 @@ psArray *flux; ///< Fluxes for possible stamps (or NULL) int footprint; ///< Half-size of stamps - psKernel *window; ///< window function generated from ensemble of stamps float normFrac; ///< Fraction of flux in window for normalisation window - int normWindow; ///< Size of window for measuring normalisation + 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 sysErr; ///< Systematic error float skyErr; ///< increase effective readnoise @@ -195,3 +197,8 @@ bool pmSubtractionStampsResetStatus (pmSubtractionStampList *stamps); + +bool pmSubtractionKernelPenaltiesStamp(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels); +bool pmSubtractionKernelPenalties(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int index); +float pmSubtractionKernelPenaltySingle(psKernel *kernel); + #endif Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c (revision 28865) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c (revision 28866) @@ -21,4 +21,5 @@ #include "pmVisual.h" +#include "pmVisualUtils.h" #include "pmHDU.h" @@ -61,5 +62,9 @@ * @return true for success */ bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) { - if (!pmVisualIsVisual() || !plotConvKernels) return true; + + if (!pmVisualTestLevel("ppsub.kernels", 1)) return true; + + if (!plotConvKernels) return true; + if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) { return false; @@ -75,5 +80,9 @@ @return true for success */ bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) { - if (!pmVisualIsVisual() || !plotStamps) return true; + + if (!pmVisualTestLevel("ppsub.stamps", 1)) return true; + + if (!plotStamps) return true; + if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) { return false; @@ -145,5 +154,8 @@ bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) { - if (!pmVisualIsVisual() || !plotLeastSquares) return true; + if (!pmVisualTestLevel("ppsub.chisq", 1)) return true; + + if (!plotLeastSquares) return true; + if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) { return false; @@ -204,5 +216,9 @@ bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) { - if (!pmVisualIsVisual() || !plotImage) return true; + + if (!pmVisualTestLevel("ppsub.images.sub", 1)) return true; + + if (!plotImage) return true; + if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) { return false; @@ -218,5 +234,6 @@ bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("ppsub.kern.final", 1)) return true; + if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) { return false; @@ -264,5 +281,6 @@ bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("ppsub.basis", 1)) return true; + if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) { return false; @@ -425,5 +443,5 @@ bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("ppsub.fit", 1)) return true; // generate 4 storage images large enough to hold the N stamps: @@ -462,5 +480,5 @@ bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("ppsub.stamp", 1)) return true; double sum; @@ -543,5 +561,6 @@ } - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("ppsub.fit", 1)) return true; + if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false; if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false; @@ -605,5 +624,6 @@ Graphdata graphdata; - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("ppsub.fit", 1)) return true; + if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;