IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29148


Ignore:
Timestamp:
Sep 12, 2010, 12:38:23 PM (16 years ago)
Author:
eugene
Message:

re-scale kernels by fwhm2 to avoid excessive dynamic range; modify residual stats to be more sensible; add threshold to min significant element of SVD inversion (1e-6); adjust penalties again (still not clear what is best...)

Location:
branches/eam_branches/ipp-20100823/psModules/src/imcombine
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtraction.c

    r29004 r29148  
    540540
    541541    psImage *conv = psImageConvolveFFT(NULL, image->image, NULL, 0, kernel); // Convolved image
     542
     543    // note: do not attempt to renormalize kernels here: cannot have different stars with
     544    // different kernel ratios
     545
    542546    int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
    543547    psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.c

    r27086 r29148  
    299299                         kernels->rms);
    300300
    301         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
    302         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
    303         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
    304         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
    305         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
    306         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
    307 
    308         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
    309         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
    310         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
    311         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
    312         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
    313         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
     301        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
     302        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
     303        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
     304        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
     305        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
     306        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
     307
     308        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
     309        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
     310        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
     311        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
     312        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
     313        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
    314314    }
    315315

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionAnalysis.h

    r26893 r29148  
    2424#define PM_SUBTRACTION_ANALYSIS_DECONV_MAX   "SUBTRACTION.DECONV.MAX"   // Maximum deconvolution fraction
    2525
    26 #define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN  "SUBTRACTION.RES.FSIGMA.MEAN"      // RMS stamp deviation
    27 #define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV "SUBTRACTION.RES.FSIGMA.STDEV"      // RMS stamp deviation
    28 #define PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN    "SUBTRACTION.RES.FMIN.MEAN"      // RMS stamp deviation
    29 #define PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV   "SUBTRACTION.RES.FMIN.STDEV"      // RMS stamp deviation
    30 #define PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN    "SUBTRACTION.RES.FMAX.MEAN"      // RMS stamp deviation
    31 #define PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV   "SUBTRACTION.RES.FMAX.STDEV"      // RMS stamp deviation
     26#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN  "SUBTRACTION.FRES.MEAN" // RMS stamp deviation
     27#define PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV "SUBTRACTION.FRES.STDEV" // RMS stamp deviation
     28#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN  "SUBTRACTION.FRES.OUTER.MEAN" // RMS stamp deviation
     29#define PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV "SUBTRACTION.FRES.OUTER.STDEV" // RMS stamp deviation
     30#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN  "SUBTRACTION.FRES.TOTAL.MEAN" // RMS stamp deviation
     31#define PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV "SUBTRACTION.FRES.TOTAL.STDEV" // RMS stamp deviation
    3232
    3333// Derive QA information about the subtraction

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c

    r29126 r29148  
    972972            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    973973            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
     974               
    974975                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
    975976                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
     977
    976978                psVectorAppend(norms, stamp->norm);
     979
    977980                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    978981                numStamps++;
    … …  
    11071110        }
    11081111
    1109 #ifdef TESTING
    1110         psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
     1112#if 0
     1113        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
     1114        psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
    11111115        psVectorWriteFile("sumVector.dat", sumVector);
     1116        psFree (save);
    11121117#endif
    11131118
    … …  
    11911196            sumVector->data.F64[normIndex] = 0.0;
    11921197
    1193             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    1194 
     1198// save the matrix and vector after the NULLs have been set
     1199#if 0
     1200            psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
     1201            psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
     1202            psVectorWriteFile("sumVector.dat", sumVector);
     1203            psFree (save);
     1204#endif
     1205
     1206            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6);
     1207            // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 3e-4);
     1208            // psVectorCopy (solution, sumVector, PS_TYPE_F64);
     1209            // psMatrixGJSolve(sumMatrix, solution);
    11951210            solution->data.F64[normIndex] = normValue;
    11961211        }
    … …  
    12651280}
    12661281
    1267 bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
    1268 
    1269     // XXX measure some useful stats on the residuals
     1282// measure some useful stats on the stamp residuals:
     1283// fResSigma : the residual stdev / total flux
     1284// fResOuter : the residual fabs / total flux for R > 2 pix
     1285// fResTotal : the residual fabs / total flux for R > 0 pix
     1286bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
     1287
    12701288    float sum = 0.0;
    12711289    float peak = 0.0;
    … …  
    12771295    }
    12781296
    1279     // only count pixels with more than X% of the source flux
    1280     // calculate stdev(dflux)
     1297    // init counters
     1298    int npix = 0;
    12811299    float dflux1 = 0.0;
    12821300    float dflux2 = 0.0;
    1283     int npix = 0;
    1284 
    1285     float dmax = 0.0;
    1286     float dmin = 0.0;
    1287 
    1288     // XXX update these with a bit more rigour
     1301    float dOuter = 0.0;
     1302    float dTotal = 0.0;
     1303
    12891304    for (int y = - footprint; y <= footprint; y++) {
    12901305        for (int x = - footprint; x <= footprint; x++) {
    1291           // float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
    1292             // if (dflux < 0.02*sum) continue;
    12931306            dflux1 += residual->kernel[y][x];
    12941307            dflux2 += PS_SQR(residual->kernel[y][x]);
    1295             // dmax = PS_MAX(residual->kernel[y][x], dmax);
    1296             // dmin = PS_MIN(residual->kernel[y][x], dmin);
    1297             dmax += fabs(residual->kernel[y][x]);
    1298 
     1308            dTotal += fabs(residual->kernel[y][x]);
    12991309            if (hypot(x,y) > 2.0) {
    1300               dmin += fabs(residual->kernel[y][x]);
     1310              dOuter += fabs(residual->kernel[y][x]);
    13011311            }
    13021312            npix ++;
    … …  
    13051315    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
    13061316    if (!isfinite(sum))  return false;
    1307     if (!isfinite(dmax)) return false;
    1308     if (!isfinite(dmin)) return false;
    13091317    if (!isfinite(peak)) return false;
    1310 
    1311     fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/sum, dmin/sum);
    1312     psVectorAppend(fSigRes, sigma/sum);
    1313     psVectorAppend(fMaxRes, dmax/sum);
    1314     psVectorAppend(fMinRes, dmin/sum);
     1318    if (!isfinite(dOuter)) return false;
     1319    if (!isfinite(dTotal)) return false;
     1320
     1321    fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum);
     1322    psVectorAppend(fResSigma, sigma/sum);
     1323    psVectorAppend(fResOuter, dOuter/sum);
     1324    psVectorAppend(fResTotal, dTotal/sum);
    13151325    return true;
    13161326}
    … …  
    13351345    pmSubtractionVisualShowFitInit (stamps);
    13361346
    1337     psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    1338     psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    1339     psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1347    psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1348    psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1349    psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    13401350
    13411351    // we want to save the residual images for the 9 brightest stamps.
    … …  
    14491459            }
    14501460
    1451             pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
     1461            pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint);
    14521462
    14531463        } else {
    … …  
    14851495            }
    14861496
    1487             pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
     1497            pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint);
    14881498        }
    14891499
    … …  
    15641574
    15651575        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    1566         psVectorStats (stats, fSigRes, NULL, NULL, 0);
    1567         kernels->fSigResMean = stats->robustMedian;
    1568         kernels->fSigResStdev = stats->robustStdev;
     1576        psVectorStats (stats, fResSigma, NULL, NULL, 0);
     1577        kernels->fResSigmaMean = stats->robustMedian;
     1578        kernels->fResSigmaStdev = stats->robustStdev;
    15691579
    15701580        psStatsInit (stats);
    1571         psVectorStats (stats, fMaxRes, NULL, NULL, 0);
    1572         kernels->fMaxResMean = stats->robustMedian;
    1573         kernels->fMaxResStdev = stats->robustStdev;
     1581        psVectorStats (stats, fResOuter, NULL, NULL, 0);
     1582        kernels->fResOuterMean = stats->robustMedian;
     1583        kernels->fResOuterStdev = stats->robustStdev;
    15741584
    15751585        psStatsInit (stats);
    1576         psVectorStats (stats, fMinRes, NULL, NULL, 0);
    1577         kernels->fMinResMean = stats->robustMedian;
    1578         kernels->fMinResStdev = stats->robustStdev;
     1586        psVectorStats (stats, fResTotal, NULL, NULL, 0);
     1587        kernels->fResTotalMean = stats->robustMedian;
     1588        kernels->fResTotalStdev = stats->robustStdev;
    15791589
    15801590        // XXX save these values somewhere
    1581         psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",
    1582                  kernels->fSigResMean, kernels->fSigResStdev,
    1583                  kernels->fMaxResMean, kernels->fMaxResStdev,
    1584                  kernels->fMinResMean, kernels->fMinResStdev);
    1585 
    1586         psFree (fSigRes);
    1587         psFree (fMaxRes);
    1588         psFree (fMinRes);
     1591        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f",
     1592                 kernels->fResSigmaMean, kernels->fResSigmaStdev,
     1593                 kernels->fResOuterMean, kernels->fResOuterStdev,
     1594                 kernels->fResTotalMean, kernels->fResTotalStdev);
     1595
     1596        psFree (fResSigma);
     1597        psFree (fResOuter);
     1598        psFree (fResTotal);
    15891599        psFree (stats);
    15901600    }
    … …  
    15921602    psFree(residual);
    15931603    psFree(polyValues);
    1594 
    15951604
    15961605    return deviations;

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c

    r29004 r29148  
    220220            // Re-normalize
    221221            // scale2D  = 1.0 / fabs(sum);
    222             scale2D  = 1.0 / sqrt(sum2);
     222            scale2D  = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
    223223            zeroNull = true;
    224224        } else {
    225225            // Odd functions: choose normalisation so that parameters have about the same strength as for even
    226226            // functions, no subtraction of null pixel because the sum is already (near) zero
    227             scale2D = 1.0 / sqrt(sum2);
     227            scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
    228228            zeroNull = false;
    229229        }
    … …  
    235235    if (forceZeroNull) {
    236236        // Force rescaling and subtraction of null pixel even though the order doesn't indicate it's even
    237         scale2D = 1.0 / fabs(sum);
     237        scale2D = 1.0 / fabs(sum) / PS_SQR(fwhm);
    238238        zeroNull = true;
    239239    }
    240240    if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
    241241        // Odd function
    242         scale2D = 1.0 / sqrt(sum2);
     242        scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
    243243    }
    244244
    … …  
    255255    if (zeroNull) {
    256256        // preCalc->kernel->kernel[0][0] -= 1.0;
    257         preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
     257        preCalc->kernel->kernel[0][0] -= sum * scale2D;
    258258    }
    259259
    … …  
    603603    kernels->sampleStamps = NULL;
    604604
    605     kernels->fSigResMean  = NAN;
    606     kernels->fSigResStdev = NAN;
    607     kernels->fMaxResMean  = NAN;
    608     kernels->fMaxResStdev = NAN;
    609     kernels->fMinResMean  = NAN;
    610     kernels->fMinResStdev = NAN;
     605    kernels->fResSigmaMean  = NAN;
     606    kernels->fResSigmaStdev = NAN;
     607    kernels->fResOuterMean  = NAN;
     608    kernels->fResOuterStdev = NAN;
     609    kernels->fResTotalMean  = NAN;
     610    kernels->fResTotalStdev = NAN;
    611611
    612612    return kernels;

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.h

    r29004 r29148  
    5151    float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
    5252    int numStamps;                      ///< Number of good stamps
    53     float fSigResMean;                  ///< mean fractional stdev of residuals
    54     float fSigResStdev;                 ///< stdev of fractional stdev of residuals
    55     float fMaxResMean;                  ///< mean fractional positive swing in residuals
    56     float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
    57     float fMinResMean;                  ///< mean fractional negative swing in residuals
    58     float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
     53    float fResSigmaMean;                ///< mean fractional stdev of residuals
     54    float fResSigmaStdev;               ///< stdev of fractional stdev of residuals
     55    float fResOuterMean;                ///< mean fractional positive swing in residuals
     56    float fResOuterStdev;               ///< stdev of fractional positive swing in residuals
     57    float fResTotalMean;                ///< mean fractional negative swing in residuals
     58    float fResTotalStdev;               ///< stdev of fractional negative swing in residuals
    5959    psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
    6060} pmSubtractionKernels;

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c

    r29126 r29148  
    919919        float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull);
    920920
    921         penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment
    922         penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment
     921        if (1) {
     922            penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment
     923            penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment
     924            // penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment
     925            // penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment
     926        } else {
     927            penalty1 = PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment
     928            penalty2 = PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment
     929        }
    923930    }
    924931    kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c

    r29126 r29148  
    152152
    153153/** Plot the least-squares matrix of each stamp */
    154 bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) {
     154bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps) {
    155155
    156156    if (!pmVisualTestLevel("ppsub.chisq", 1)) return true;
    … …  
    209209    pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
    210210
     211    if (0) {
     212        static int count = 0;
     213        char filename[64];
     214        sprintf (filename, "chisq.%02d.fits", count);
     215        count ++;
     216        psFits *fits = psFitsOpen (filename, "w");
     217        psFitsWriteImage (fits, NULL, canvas32, 0, NULL);
     218        psFitsClose (fits);
     219    }
     220
    211221    pmVisualAskUser(&plotLeastSquares);
    212222    psFree(canvas);
    … …  
    299309        if (!isfinite(stamp->flux)) continue;
    300310        if (!stamp->convolutions1 && !stamp->convolutions2) continue;
     311        fprintf (stderr, "flux: %f, maxFlux: %f  ", stamp->flux, maxFlux);
    301312        if (!maxStamp) {
    302313            maxFlux = stamp->flux;
    303314            maxStamp = stamp;
     315            fprintf (stderr, "maxStamp %d\n", i);
    304316            continue;
     317        } else {
     318            fprintf (stderr, "\n");
    305319        }
    306320        if (stamp->flux > maxFlux) {
    … …  
    339353           
    340354            double sum = 0.0;
     355            double sum2 = 0.0;
    341356            for (int y = -footprint; y <= footprint; y++) {
    342357                for (int x = -footprint; x <= footprint; x++) {
    343358                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
    344359                    sum += kernel->kernel[y][x];
     360                    sum2 += PS_SQR(kernel->kernel[y][x]);
    345361                }
    346362            }
    347             fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     363            fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
    348364        }               
    349365        pmVisualScaleImage(kapa2, output, "Image", 0, true);
    350     }                                   
     366
     367        if (0) {
     368            psFits *fits = psFitsOpen("basis.1.fits", "w");
     369            psFitsWriteImage(fits, NULL, output, 0, NULL);
     370            psFitsClose(fits);
     371        }
     372    }
    351373       
    352374    if (maxStamp->convolutions2) {
    … …  
    373395           
    374396            double sum = 0.0;
     397            double sum2 = 0.0;
    375398            for (int y = -footprint; y <= footprint; y++) {
    376399                for (int x = -footprint; x <= footprint; x++) {
    377400                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
    378401                    sum += kernel->kernel[y][x];
     402                    sum2 += PS_SQR(kernel->kernel[y][x]);
    379403                }
    380404            }
    381             fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     405            fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
    382406        }               
    383407        pmVisualScaleImage(kapa2, output, "Image", 1, true);
     408
     409        if (0) {
     410            psFits *fits = psFitsOpen("basis.2.fits", "w");
     411            psFitsWriteImage(fits, NULL, output, 0, NULL);
     412            psFitsClose(fits);
     413        }
    384414    }                                   
    385415       
Note: See TracChangeset for help on using the changeset viewer.