Changeset 30266

Timestamp:

Jan 12, 2011, 9:32:08 PM (15 years ago)

Author:

eugene

Message:

add errors to kernel coeffs; calculate and plot chisq and moments

Location:

branches/eam_branches/ipp-20101205/psModules/src

Files:

• extras/pmVisual.c (modified) (1 diff)
• extras/pmVisual.h (modified) (1 diff)
• imcombine/pmSubtractionEquation.c (modified) (18 diffs)
• imcombine/pmSubtractionEquation.h (modified) (1 diff)
• imcombine/pmSubtractionKernels.c (modified) (3 diffs)
• imcombine/pmSubtractionKernels.h (modified) (1 diff)
• imcombine/pmSubtractionMatch.c (modified) (1 diff)
• imcombine/pmSubtractionStamps.c (modified) (1 diff)
• imcombine/pmSubtractionVisual.c (modified) (2 diffs)
• imcombine/pmSubtractionVisual.h (modified) (2 diffs)

branches/eam_branches/ipp-20101205/psModules/src/extras/pmVisual.c

@@ -306 +306 @@
 }
 
+bool pmVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max) {
+
+    KiiImage image;
+    KapaImageData data;
+    Coords coords;
+
+    strcpy (coords.ctype, "RA---TAN");
+
+    image.data2d = inImage->data.F32;
+    image.Nx = inImage->numCols;
+    image.Ny = inImage->numRows;
+
+    strcpy (data.name, name);
+    strcpy (data.file, name);
+    data.zero = min;
+    data.range = max - min;
+    data.logflux = 0;
+
+    KiiSetChannel (kapaFD, channel);
+    KiiNewPicture2D (kapaFD, &image, &data, &coords);
+
+    return true;
+}
 
 psImage* pmVisualImageToFloat(psImage *image) {

branches/eam_branches/ipp-20101205/psModules/src/extras/pmVisual.h

@@ -64 +64 @@
                         const char *name, int channel, bool clip);
 
+bool pmVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
 
 /** Calculate statistics on an image.

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c

@@ -8 +8 @@
 
 #include "pmErrorCodes.h"
+#include "pmVisual.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
@@ -16 +17 @@
 #include "pmSubtractionVisual.h"
 
-//#define TESTING                         // TESTING output for debugging; may not work with threads!
-
-//#define USE_WEIGHT                      // Include weight (1/variance) in equation?
-
-// XXX TEST:
-//# define USE_WINDOW                      // window to avoid neighbor contamination
+//# define TESTING                         // TESTING output for debugging; may not work with threads!
+//# define USE_WEIGHT                      // Include weight (1/variance) in equation?
+# define USE_WINDOW                      // window to avoid neighbor contamination
+// XXX if we want to have a weight and window, we'll need to pass through to the chisq function
 
 # define PENALTY false
@@ -983 +982 @@
     }
 
-    pmSubtractionVisualPlotLeastSquares(stamps);
     pmSubtractionVisualShowKernels((pmSubtractionKernels  *)kernels);
     pmSubtractionVisualShowBasis(stamps);
+    pmSubtractionVisualPlotLeastSquares(stamps);
 
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
@@ -1080 +1079 @@
         }
 
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
+
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
@@ -1087 +1088 @@
 #endif
 
+        psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
         solution = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -1096 +1098 @@
         if (1) {
             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+            invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
         } else {
             psVector *PERM = NULL;
@@ -1104 +1107 @@
         }
 
-# if (0)
+# if (1)
         for (int i = 0; i < solution->n; i++) {
-            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]);
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
         }
 # endif
@@ -1113 +1116 @@
             kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
             psVectorInit(kernels->solution1, 0.0);
+            kernels->solution1err = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+            psVectorInit(kernels->solution1err, 0.0);
         }
 
@@ -1120 +1125 @@
         for (int i = 0; i < numKernels * numPoly; i++) {
             kernels->solution1->data.F64[i] = solution->data.F64[i];
+            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]);
         }
 
@@ -1131 +1137 @@
         psFree(sumVector);
         psFree(sumMatrix);
+        psFree(invMatrix);
 
     } else {
@@ -1160 +1167 @@
         }
 
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
+
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
@@ -1171 +1180 @@
         }
 
+        psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
         solution = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -1176 +1186 @@
 
         // DUAL solution
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-
-#if (0)
+# if (1)
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-7);
+        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
+# endif
+# if (0)
+        psMatrixLUSolve(sumMatrix, sumVector);
+        solution = sumVector;
+        invMatrix = sumMatrix;
+# endif
+
+#if (1)
         for (int i = 0; i < solution->n; i++) {
-            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
+            fprintf(stderr, "Dual solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
         }
 #endif
@@ -1187 +1205 @@
             kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64);
             psVectorInit (kernels->solution1, 0.0);
+            kernels->solution1err = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+            psVectorInit(kernels->solution1err, 0.0);
         }
         if (!kernels->solution2) {
             kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
             psVectorInit (kernels->solution2, 0.0);
+            kernels->solution2err = psVectorAlloc(numSolution2, PS_TYPE_F64);
+            psVectorInit(kernels->solution2err, 0.0);
         }
 
@@ -1205 +1227 @@
             kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue;
             kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
+
+            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]) * stamps->normValue;
+            int i2 = i + numSolution1;
+            kernels->solution2err->data.F64[i] = sqrt(invMatrix->data.F64[i2][i2]);
         }
 
@@ -1213 +1239 @@
         kernels->solution1->data.F64[bgIndex] = 0.0;
 
+        psFree(solution);
+        psFree(sumVector);
         psFree(sumMatrix);
-        psFree(sumVector);
-        psFree(solution);
+        psFree(invMatrix);
     }
 
@@ -1224 +1251 @@
     if (psTraceGetLevel("psModules.imcombine") >= 7) {
         for (int i = 0; i < kernels->solution1->n; i++) {
-            psTrace("psModules.imcombine", 7, "Solution 1 %d: %f\n", i, kernels->solution1->data.F64[i]);
+            psTrace("psModules.imcombine", 7, "Solution 1 %d: %f +/- %f\n", i, kernels->solution1->data.F64[i], kernels->solution1err->data.F64[i]);
         }
         if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
             for (int i = 0; i < kernels->solution2->n; i++) {
-                psTrace("psModules.imcombine", 7, "Solution 2 %d: %f\n", i, kernels->solution2->data.F64[i]);
+                psTrace("psModules.imcombine", 7, "Solution 2 %d: %f +/- %f\n", i, kernels->solution2->data.F64[i], kernels->solution2err->data.F64[i]);
             }
         }
@@ -1280 +1307 @@
     psVectorAppend(fResOuter, dOuter/sum);
     psVectorAppend(fResTotal, dTotal/sum);
+    return true;
+}
+
+// given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) {
+
+    int npix = 0;
+    float chisq = 0;
+
+    // get the chisq
+    for (int y = residual->yMin; y <= residual->yMax; y++) {
+        for (int x = residual->xMin; x <= residual->xMax; x++) {
+            chisq  += PS_SQR(residual->kernel[y][x]);
+            npix ++;
+        }
+    }
+    psVectorAppend(chisqVector, chisq / npix);
+
+    float value1 = 0;
+    float value2 = 0;
+    float flux2 = 0;
+    float fluxX = 0;
+    float fluxY = 0;
+    float fluxX2 = 0;
+    float fluxY2 = 0;
+
+    float fluxC1 = 0;
+    float fluxC2 = 0;
+
+    float moment = 0;
+
+    // get the moments from convolved1
+    if (convolved1) {
+        for (int y = residual->yMin; y <= residual->yMax; y++) {
+            for (int x = residual->xMin; x <= residual->xMax; x++) {
+                value1  = convolved1->kernel[y][x];
+                value2  = PS_SQR(value1);
+
+                if (weight) {
+                    value2 *= weight->kernel[y][x];
+                }
+                if (window) {
+                    value1 *= window->kernel[y][x];
+                    value2 *= window->kernel[y][x];
+                }
+
+                fluxC1 += value1;
+                flux2  += value2;
+                fluxX  += x * value2;
+                fluxY  += y * value2;
+                fluxX2 += PS_SQR(x) * value2;
+                fluxY2 += PS_SQR(y) * value2;
+            }
+        }
+        // float Mx = fluxX / flux2;
+        // float My = fluxY / flux2;
+        float Mxx = fluxX2 / flux2;
+        float Myy = fluxY2 / flux2;
+
+        // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
+        moment += Mxx + Myy;
+    }
+
+    // get the moments from convolved1
+    if (convolved2) {
+        for (int y = residual->yMin; y <= residual->yMax; y++) {
+            for (int x = residual->xMin; x <= residual->xMax; x++) {
+                value1  = convolved2->kernel[y][x];
+                value2  = PS_SQR(value1);
+
+                if (weight) {
+                    value2 *= weight->kernel[y][x];
+                }
+                if (window) {
+                    value1 *= window->kernel[y][x];
+                    value2 *= window->kernel[y][x];
+                }
+
+                fluxC2 += value1;
+                flux2  += value2;
+                fluxX  += x * value2;
+                fluxY  += y * value2;
+                fluxX2 += PS_SQR(x) * value2;
+                fluxY2 += PS_SQR(y) * value2;
+            }
+        }
+        // float Mx = fluxX / flux2;
+        // float My = fluxY / flux2;
+        float Mxx = fluxX2 / flux2;
+        float Myy = fluxY2 / flux2;
+
+        // fprintf (stderr, "conv2, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
+        moment += Mxx + Myy;
+    }
+
+    float flux = fluxC1 + fluxC2;
+    
+    if (convolved1 && convolved2) {
+        moment *= 0.5;
+        flux *= 0.5;
+    }
+    psVectorAppend(momentVector, moment);
+    psVectorAppend(fluxesVector, flux);
+    return true;
+}
+
+// typedef struct {
+//     float moment;
+//     float chisq;
+// } diffStat;
+
+bool pmSubtractionCalculateChisqAndMoments(pmSubtractionStampList *stamps,
+                                           pmSubtractionKernels *kernels)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
+
+    psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *chisq = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+
+    int footprint = stamps->footprint; // Half-size of stamps
+    int numKernels = kernels->num;      // Number of kernels
+
+    psImage *polyValues = NULL;         // Polynomial values
+
+    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+
+    // storage for the convolved source image (only convolved1 is used in SINGLE mode)
+    psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            continue;
+        }
+
+        // Calculate coefficients of the kernel basis functions
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+
+        // Calculate residuals
+        psImageInit(residual->image, 0.0);
+
+    psKernel *weight = NULL;
+    psKernel *window = NULL;
+
+#ifdef USE_WEIGHT
+    weight = stamp->weight;
+#endif
+#ifdef USE_WINDOW
+    window = stamps->window;
+#endif
+
+        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
+
+            // the single-direction psf match code attempts to find the kernel such that:
+            // source * kernel = target.  we need to assign 'source' and 'target' correctly
+            // depending on which of image1 or image2 we asked to be convolved.
+
+            psKernel *target;           // Target postage stamp (convolve source to match the target)
+            psKernel *source;           // Source postage stamp (convolve source to match the target)
+            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
+
+            // init the accumulation image
+            psImageInit(convolved1->image, 0.0);
+
+            switch (kernels->mode) {
+              case PM_SUBTRACTION_MODE_1:
+                target = stamp->image2;
+                source = stamp->image1;
+                convolutions = stamp->convolutions1;
+                break;
+              case PM_SUBTRACTION_MODE_2:
+                target = stamp->image1;
+                source = stamp->image2;
+                convolutions = stamp->convolutions2;
+                break;
+              default:
+                psAbort("Unsupported subtraction mode: %x", kernels->mode);
+            }
+
+            // generate the convolved source image (sum over kernels)
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *convolution = convolutions->data[j]; // Convolution
+                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        convolved1->kernel[y][x] += convolution->kernel[y][x] * coefficient;
+                    }
+                }
+            }
+
+            // generate the residual image
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    convolved1->kernel[y][x] += source->kernel[y][x] * norm;
+                    residual->kernel[y][x] = target->kernel[y][x] - convolved1->kernel[y][x] - background;
+                }
+            }
+            // XXX if we want to have a weight and window, we'll need to pass through to here
+            pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, NULL, residual, weight, window);
+
+        } else {
+
+            // Dual convolution
+            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
+            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
+            psKernel *image1 = stamp->image1; // The first image
+            psKernel *image2 = stamp->image2; // The second image
+
+            // init the accumulation images
+            psImageInit(convolved1->image, 0.0);
+            psImageInit(convolved2->image, 0.0);
+
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
+                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
+                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
+                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
+
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        // NOTE sign for coeff2
+                        convolved1->kernel[y][x] += +conv1->kernel[y][x] * coeff1;
+                        convolved2->kernel[y][x] += -conv2->kernel[y][x] * coeff2;
+                    }
+                }
+            }
+
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    convolved1->kernel[y][x] += image1->kernel[y][x] * norm;
+                    convolved2->kernel[y][x] += image2->kernel[y][x];
+                    residual->kernel[y][x] = convolved2->kernel[y][x] - convolved1->kernel[y][x] - background;
+                }
+            }
+
+            if (0) {
+                psFitsWriteImageSimple("conv1.fits", convolved1->image, NULL);
+                psFitsWriteImageSimple("conv2.fits", convolved2->image, NULL);
+                psFitsWriteImageSimple("resid.fits", residual->image,   NULL);
+                pmVisualAskUser(NULL);
+            } 
+
+            pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, convolved2, residual, weight, window);
+        }
+    }
+
+    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments);
+
+    // find the mean chisq and mean moment
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+    psVectorStats (stats, chisq, NULL, NULL, 0);
+    float chisqValue = stats->sampleMean;
+
+    psStatsInit(stats);
+    psVectorStats (stats, moments, NULL, NULL, 0);
+    float momentValue = stats->sampleMean;
+
+    fprintf (stderr, "chisq: %f, moment: %f\n", chisqValue, momentValue);
+
+    psFree(stats);
+    psFree(chisq);
+    psFree(fluxes);
+    psFree(moments);
+
+    psFree(residual);
+    psFree(convolved1);
+    psFree(convolved2);
+    psFree(polyValues);
+
     return true;
 }

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h

@@ -92 +92 @@
 bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window);
 
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window);
+
+bool pmSubtractionCalculateChisqAndMoments(pmSubtractionStampList *stamps,
+                                           pmSubtractionKernels *kernels);
+
 #endif

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.c

@@ -30 +30 @@
     psFree(kernels->solution1);
     psFree(kernels->solution2);
+    psFree(kernels->solution1err);
+    psFree(kernels->solution2err);
     psFree(kernels->sampleStamps);
 }
@@ -745 +747 @@
     kernels->solution1 = NULL;
     kernels->solution2 = NULL;
+    kernels->solution1err = NULL;
+    kernels->solution2err = NULL;
     kernels->mean = NAN;
     kernels->rms = NAN;
@@ -1435 +1439 @@
     out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
     out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
+    out->solution1err = in->solution1err ? psVectorCopy(NULL, in->solution1err, PS_TYPE_F64) : NULL;
+    out->solution2err = in->solution2err ? psVectorCopy(NULL, in->solution2err, PS_TYPE_F64) : NULL;
     out->sampleStamps = psMemIncrRefCounter(in->sampleStamps);
 

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionKernels.h

@@ -48 +48 @@
     pmSubtractionMode mode;             ///< Mode for subtraction
     psVector *solution1, *solution2;    ///< Solution for the PSF matching
+    psVector *solution1err, *solution2err; ///< error in solution for the PSF matching
     // Quality information
     float mean, rms;                    ///< Mean and RMS of chi^2 from stamps

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c

@@ -743 +743 @@
                 memCheck("   calculate deviations");
 
+                pmSubtractionCalculateChisqAndMoments(stamps, kernels); // Stamp deviations
+
                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
                 numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c

@@ -812 +812 @@
     stamps->normWindow2 = 2.0*R2;
     psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
+
+    // Generate a weighting window based on the kron radii
+    { 
+        float radius = 1.5 * PS_MAX(R1, R2);
+
+        psImageInit(stamps->window->image, 0.0);
+
+        for (int y = -size; y <= size; y++) {
+            for (int x = -size; x <= size; x++) {
+                if (hypot(x,y) > radius) continue;
+                stamps->window->kernel[y][x] = 1.0;
+            }
+        }
+    }
 
 # else

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c

@@ -210 +210 @@
         psImageOverlaySection(canvas, im, x0, y0, "=");
         stampNum++;
+
+        // renormalize the section
+        float maxValue = 0;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]);
+            }
+        }
+        if (maxValue == 0.0) continue;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                canvas->data.F64[y0 + iy][x0 + ix] /= maxValue;
+            }
+        }
     }
 
     psImage *canvas32 = pmVisualImageToFloat(canvas);
     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);
+    psFree(canvas32);
+    return true;
+}
+
+/** Plot the least-squares matrix of each stamp */
+bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed) {
+
+    if (!pmVisualTestLevel("ppsub.chisq", 1)) return true;
+
+    if (!plotLeastSquares) return true;
+
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
+        return false;
+    }
+
+    psImage *matrixNorm = psImageCopy(NULL, matrixIn, PS_TYPE_F64);
+    {
+        // renormalize the matrix
+        float maxValue = 0;
+        for (int iy = 0; iy < matrixNorm->numRows; iy++) {
+            for (int ix = 0; ix < matrixNorm->numCols; ix++) {
+                maxValue = PS_MAX(maxValue, matrixNorm->data.F64[iy][ix]);
+            }
+        }
+        for (int iy = 0; iy < matrixNorm->numRows; iy++) {
+            for (int ix = 0; ix < matrixNorm->numCols; ix++) {
+                matrixNorm->data.F64[iy][ix] /= maxValue;
+            }
+        }
+    }
+
+    // Find the stamp size
+    int imageMax = -1;
+    int numStamps = 0;
+    psElemType type = PS_TYPE_F64;
+
+    for(int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i];
+        if (stamp == NULL) continue;
+
+        psImage *im = stamp->matrix;
+        if (im == NULL) continue;
+
+        imageMax = PS_MAX(imageMax, im->numCols);
+        imageMax = PS_MAX(imageMax, im->numRows);
+        numStamps++;
+        type = im->type.type;
+    }
+    if (imageMax == -1) return false;
+
+    int border = 15;
+    imageMax += border;
+    int tileRowCount = (int) ceil(sqrt(numStamps));
+    int canvasX = tileRowCount * (imageMax);
+    int canvasY = tileRowCount * (imageMax);
+    psImage *canvas = psImageAlloc (canvasX, canvasY, type);
+    psImageInit (canvas, NAN);
+
+    // overlay the images
+    int stampNum = 0;
+    int stampListNum = 0;
+    while (stampNum < numStamps) {
+        int x0 = (imageMax) * (stampNum % tileRowCount);
+        int y0 = (imageMax) * (stampNum / tileRowCount);
+
+        pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++];
+        if (stamp == NULL) continue;
+
+        psImage *im = stamp->matrix;
+        if (im == NULL) continue;
+
+        stampNum++;
+
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+
+        // renormalize the section
+        float maxValue = 0;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]);
+            }
+        }
+        if (maxValue == 0.0) continue;
+        for (int iy = 0; iy < im->numRows; iy++) {
+            for (int ix = 0; ix < im->numCols; ix++) {
+                canvas->data.F64[y0 + iy][x0 + ix] = (im->data.F64[iy][ix] / maxValue - matrixNorm->data.F64[iy][ix]) / matrixNorm->data.F64[iy][ix];
+            }
+        }
+    }
+
+    psImage *canvas32 = pmVisualImageToFloat(canvas);
+    pmVisualRangeImage(kapa2, canvas32, "Least_Squares", 0, -100.0, 100.0);
 
     if (0) {
@@ -722 +841 @@
 }
 
+// plot log(flux) vs log(chisq), log(flux) vs log(moments), log(chisq) vs log(moments)
+bool pmSubtractionVisualPlotChisqAndMoments(psVector *fluxes, psVector *chisq, psVector *moments) {
+
+    Graphdata graphdata;
+    KapaSection section;
+
+    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
+
+    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;
+
+    KapaClearSections (kapa3);
+    KapaInitGraph (&graphdata);
+    KiiResize(kapa3, 1500, 500);
+
+    psVector *lchi = psVectorAlloc (fluxes->n, PS_TYPE_F32);
+    psVector *lflx = psVectorAlloc (fluxes->n, PS_TYPE_F32);
+    psVector *lMxx = psVectorAlloc (fluxes->n, PS_TYPE_F32);
+
+    // construct the plot vectors
+    for (int i = 0; i < fluxes->n; i++) {
+        lchi->data.F32[i] = log10(chisq->data.F32[i]);
+        lflx->data.F32[i] = log10(fluxes->data.F32[i]);
+        lMxx->data.F32[i] = log10(moments->data.F32[i]);
+    }
+
+    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
+
+    // section 1: lflux vs lchi
+    section.dx = 0.33;
+    section.dy = 1.00;
+    section.x  = 0.00;
+    section.y  = 0.00;
+    section.name = psStringCopy ("flux.v.chi");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualScaleGraphdata(&graphdata, lflx, lchi, false);
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
+    KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x");
+    KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "y");
+
+    // section 2: lflux vs lMxx
+    section.dx = 0.33;
+    section.dy = 1.00;
+    section.x  = 0.33;
+    section.y  = 0.00;
+    section.name = psStringCopy ("flux.v.mom");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualScaleGraphdata(&graphdata, lflx, lMxx, false);
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
+    KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x");
+    KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y");
+
+    // section 1: lflux vs lchi
+    section.dx = 0.33;
+    section.dy = 1.00;
+    section.x  = 0.66;
+    section.y  = 0.00;
+    section.name = psStringCopy ("chi.v.mom");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualScaleGraphdata(&graphdata, lchi, lMxx, false);
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
+    KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "x");
+    KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y");
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
 #else
 bool pmSubtractionVisualClose(void) {return true;}

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.h

@@ -6 +6 @@
 bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro);
 bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps);
+bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed);
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub);
 bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps);
@@ -13 +14 @@
 bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels);
 bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps);
+bool pmSubtractionVisualPlotChisqAndMoments(psVector *fluxes, psVector *chisq, psVector *moments);
 
 #endif