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

Timestamp:

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

Author:

eugene

Message:

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

File:

• branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c (modified) (2 diffs)

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;}