Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtraction.c (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtraction.c (revision 29148) @@ -540,4 +540,8 @@ psImage *conv = psImageConvolveFFT(NULL, image->image, NULL, 0, kernel); // Convolved image + + // note: do not attempt to renormalize kernels here: cannot have different stars with + // different kernel ratios + int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 29148) @@ -299,17 +299,17 @@ kernels->rms); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean); - psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, 0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, 0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, 0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); + + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN, 0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN, 0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN, 0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev); } Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.h =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 29148) @@ -24,10 +24,10 @@ #define PM_SUBTRACTION_ANALYSIS_DECONV_MAX "SUBTRACTION.DECONV.MAX" // Maximum deconvolution fraction -#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN "SUBTRACTION.RES.FSIGMA.MEAN" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV "SUBTRACTION.RES.FSIGMA.STDEV" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN "SUBTRACTION.RES.FMIN.MEAN" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV "SUBTRACTION.RES.FMIN.STDEV" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN "SUBTRACTION.RES.FMAX.MEAN" // RMS stamp deviation -#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV "SUBTRACTION.RES.FMAX.STDEV" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN "SUBTRACTION.FRES.MEAN" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV "SUBTRACTION.FRES.STDEV" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN "SUBTRACTION.FRES.OUTER.MEAN" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV "SUBTRACTION.FRES.OUTER.STDEV" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN "SUBTRACTION.FRES.TOTAL.MEAN" // RMS stamp deviation +#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV "SUBTRACTION.FRES.TOTAL.STDEV" // RMS stamp deviation // Derive QA information about the subtraction Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c (revision 29148) @@ -972,7 +972,10 @@ pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + psVectorAppend(norms, stamp->norm); + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; @@ -1107,7 +1110,9 @@ } -#ifdef TESTING - psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); +#if 0 + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); #endif @@ -1191,6 +1196,16 @@ sumVector->data.F64[normIndex] = 0.0; - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - +// 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; } @@ -1265,7 +1280,10 @@ } -bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { - - // XXX measure some useful stats on the residuals +// measure some useful stats on the stamp residuals: +// fResSigma : the residual stdev / total flux +// fResOuter : the residual fabs / total flux for R > 2 pix +// fResTotal : the residual fabs / total flux for R > 0 pix +bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { + float sum = 0.0; float peak = 0.0; @@ -1277,26 +1295,18 @@ } - // only count pixels with more than X% of the source flux - // calculate stdev(dflux) + // init counters + int npix = 0; float dflux1 = 0.0; float dflux2 = 0.0; - int npix = 0; - - float dmax = 0.0; - float dmin = 0.0; - - // XXX update these with a bit more rigour + float dOuter = 0.0; + float dTotal = 0.0; + for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - // float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm); - // if (dflux < 0.02*sum) continue; dflux1 += residual->kernel[y][x]; dflux2 += PS_SQR(residual->kernel[y][x]); - // dmax = PS_MAX(residual->kernel[y][x], dmax); - // dmin = PS_MIN(residual->kernel[y][x], dmin); - dmax += fabs(residual->kernel[y][x]); - + dTotal += fabs(residual->kernel[y][x]); if (hypot(x,y) > 2.0) { - dmin += fabs(residual->kernel[y][x]); + dOuter += fabs(residual->kernel[y][x]); } npix ++; @@ -1305,12 +1315,12 @@ float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix)); if (!isfinite(sum)) return false; - if (!isfinite(dmax)) return false; - if (!isfinite(dmin)) return false; if (!isfinite(peak)) return false; - - fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/sum, dmin/sum); - psVectorAppend(fSigRes, sigma/sum); - psVectorAppend(fMaxRes, dmax/sum); - psVectorAppend(fMinRes, dmin/sum); + if (!isfinite(dOuter)) return false; + if (!isfinite(dTotal)) return false; + + fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum); + psVectorAppend(fResSigma, sigma/sum); + psVectorAppend(fResOuter, dOuter/sum); + psVectorAppend(fResTotal, dTotal/sum); return true; } @@ -1335,7 +1345,7 @@ pmSubtractionVisualShowFitInit (stamps); - psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); - psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); - psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); // we want to save the residual images for the 9 brightest stamps. @@ -1449,5 +1459,5 @@ } - pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint); + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint); } else { @@ -1485,5 +1495,5 @@ } - pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint); + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint); } @@ -1564,27 +1574,27 @@ psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psVectorStats (stats, fSigRes, NULL, NULL, 0); - kernels->fSigResMean = stats->robustMedian; - kernels->fSigResStdev = stats->robustStdev; + psVectorStats (stats, fResSigma, NULL, NULL, 0); + kernels->fResSigmaMean = stats->robustMedian; + kernels->fResSigmaStdev = stats->robustStdev; psStatsInit (stats); - psVectorStats (stats, fMaxRes, NULL, NULL, 0); - kernels->fMaxResMean = stats->robustMedian; - kernels->fMaxResStdev = stats->robustStdev; + psVectorStats (stats, fResOuter, NULL, NULL, 0); + kernels->fResOuterMean = stats->robustMedian; + kernels->fResOuterStdev = stats->robustStdev; psStatsInit (stats); - psVectorStats (stats, fMinRes, NULL, NULL, 0); - kernels->fMinResMean = stats->robustMedian; - kernels->fMinResStdev = stats->robustStdev; + psVectorStats (stats, fResTotal, NULL, NULL, 0); + kernels->fResTotalMean = stats->robustMedian; + kernels->fResTotalStdev = stats->robustStdev; // XXX save these values somewhere - psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f", - kernels->fSigResMean, kernels->fSigResStdev, - kernels->fMaxResMean, kernels->fMaxResStdev, - kernels->fMinResMean, kernels->fMinResStdev); - - psFree (fSigRes); - psFree (fMaxRes); - psFree (fMinRes); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f", + kernels->fResSigmaMean, kernels->fResSigmaStdev, + kernels->fResOuterMean, kernels->fResOuterStdev, + kernels->fResTotalMean, kernels->fResTotalStdev); + + psFree (fResSigma); + psFree (fResOuter); + psFree (fResTotal); psFree (stats); } @@ -1592,5 +1602,4 @@ psFree(residual); psFree(polyValues); - return deviations; Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c (revision 29148) @@ -220,10 +220,10 @@ // Re-normalize // scale2D = 1.0 / fabs(sum); - scale2D = 1.0 / sqrt(sum2); + scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm); zeroNull = true; } else { // Odd functions: choose normalisation so that parameters have about the same strength as for even // functions, no subtraction of null pixel because the sum is already (near) zero - scale2D = 1.0 / sqrt(sum2); + scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm); zeroNull = false; } @@ -235,10 +235,10 @@ if (forceZeroNull) { // Force rescaling and subtraction of null pixel even though the order doesn't indicate it's even - scale2D = 1.0 / fabs(sum); + scale2D = 1.0 / fabs(sum) / PS_SQR(fwhm); zeroNull = true; } if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) { // Odd function - scale2D = 1.0 / sqrt(sum2); + scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm); } @@ -255,5 +255,5 @@ if (zeroNull) { // preCalc->kernel->kernel[0][0] -= 1.0; - preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2); + preCalc->kernel->kernel[0][0] -= sum * scale2D; } @@ -603,10 +603,10 @@ kernels->sampleStamps = NULL; - kernels->fSigResMean = NAN; - kernels->fSigResStdev = NAN; - kernels->fMaxResMean = NAN; - kernels->fMaxResStdev = NAN; - kernels->fMinResMean = NAN; - kernels->fMinResStdev = NAN; + kernels->fResSigmaMean = NAN; + kernels->fResSigmaStdev = NAN; + kernels->fResOuterMean = NAN; + kernels->fResOuterStdev = NAN; + kernels->fResTotalMean = NAN; + kernels->fResTotalStdev = NAN; return kernels; Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.h (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.h (revision 29148) @@ -51,10 +51,10 @@ float mean, rms; ///< Mean and RMS of chi^2 from stamps int numStamps; ///< Number of good stamps - float fSigResMean; ///< mean fractional stdev of residuals - float fSigResStdev; ///< stdev of fractional stdev of residuals - float fMaxResMean; ///< mean fractional positive swing in residuals - float fMaxResStdev; ///< stdev of fractional positive swing in residuals - float fMinResMean; ///< mean fractional negative swing in residuals - float fMinResStdev; ///< stdev of fractional negative swing in residuals + float fResSigmaMean; ///< mean fractional stdev of residuals + float fResSigmaStdev; ///< stdev of fractional stdev of residuals + float fResOuterMean; ///< mean fractional positive swing in residuals + float fResOuterStdev; ///< stdev of fractional positive swing in residuals + float fResTotalMean; ///< mean fractional negative swing in residuals + float fResTotalStdev; ///< stdev of fractional negative swing in residuals psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations } pmSubtractionKernels; Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c (revision 29148) @@ -919,6 +919,13 @@ float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull); - 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 + if (1) { + 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 + // 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 + } else { + penalty1 = PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment + penalty2 = PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment + } } kernels->penalties1->data.F32[index] = kernels->penalty * penalty1; Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c (revision 29147) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c (revision 29148) @@ -152,5 +152,5 @@ /** Plot the least-squares matrix of each stamp */ -bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) { +bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps) { if (!pmVisualTestLevel("ppsub.chisq", 1)) return true; @@ -209,4 +209,14 @@ pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true); + if (0) { + static int count = 0; + char filename[64]; + sprintf (filename, "chisq.%02d.fits", count); + count ++; + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, canvas32, 0, NULL); + psFitsClose (fits); + } + pmVisualAskUser(&plotLeastSquares); psFree(canvas); @@ -299,8 +309,12 @@ if (!isfinite(stamp->flux)) continue; if (!stamp->convolutions1 && !stamp->convolutions2) continue; + fprintf (stderr, "flux: %f, maxFlux: %f ", stamp->flux, maxFlux); if (!maxStamp) { maxFlux = stamp->flux; maxStamp = stamp; + fprintf (stderr, "maxStamp %d\n", i); continue; + } else { + fprintf (stderr, "\n"); } if (stamp->flux > maxFlux) { @@ -339,14 +353,22 @@ double sum = 0.0; + double sum2 = 0.0; for (int y = -footprint; y <= footprint; y++) { for (int x = -footprint; x <= footprint; x++) { output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; sum += kernel->kernel[y][x]; + sum2 += PS_SQR(kernel->kernel[y][x]); } } - fprintf (stderr, "kernel %d, sum %f\n", i, sum); + fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); } pmVisualScaleImage(kapa2, output, "Image", 0, true); - } + + if (0) { + psFits *fits = psFitsOpen("basis.1.fits", "w"); + psFitsWriteImage(fits, NULL, output, 0, NULL); + psFitsClose(fits); + } + } if (maxStamp->convolutions2) { @@ -373,13 +395,21 @@ double sum = 0.0; + double sum2 = 0.0; for (int y = -footprint; y <= footprint; y++) { for (int x = -footprint; x <= footprint; x++) { output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; sum += kernel->kernel[y][x]; + sum2 += PS_SQR(kernel->kernel[y][x]); } } - fprintf (stderr, "kernel %d, sum %f\n", i, sum); + fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); } pmVisualScaleImage(kapa2, output, "Image", 1, true); + + if (0) { + psFits *fits = psFitsOpen("basis.2.fits", "w"); + psFitsWriteImage(fits, NULL, output, 0, NULL); + psFitsClose(fits); + } }