Changeset 25022 for branches/eam_branches/20090715/psModules/src/imcombine

Timestamp:

Aug 7, 2009, 12:45:24 PM (17 years ago)

Author:

eugene

Message:

update from mainline

Location:

branches/eam_branches/20090715

Files:

• . (modified) (1 prop)
• psModules/src/imcombine/pmReadoutCombine.c (modified) (3 diffs)
• psModules/src/imcombine/pmReadoutCombine.h (modified) (1 diff)
• psModules/src/imcombine/pmSubtractionEquation.c (modified) (5 diffs)

branches/eam_branches/20090715/psModules/src/imcombine/pmReadoutCombine.c

@@ -281 +281 @@
     #endif
 
+    // set up windows for visualization (if selected)
+    pmReadoutCombineVisualInit();
+
     // Dereference output products
     psF32 **outputImage  = output->image->data.F32; // Output image
@@ -397 +400 @@
             outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
 
+            if (!isfinite(outputImage[yOut][xOut])) {
+                pmReadoutCombineVisualPixels(pixels, mask, outputImage[yOut][xOut]);
+            }
+
             if (isnan(outputImage[yOut][xOut])) {
                 outputImage[yOut][xOut] = NAN;
@@ -434 +441 @@
 }
 
+#if (HAVE_KAPA)
+#include <kapa.h>
+#include "pmKapaPlots.h"
+#include "pmVisual.h"
+
+static int kapa = -1;
+static bool plotFlag = true;
+
+// this init function only gets the ordinates for the first readout...
+bool pmReadoutCombineVisualInit(void) {
+    
+    if (!pmVisualIsVisual()) return true;
+
+    // skip if we have already opened the windows (or if none are requested...)
+    if (kapa != -1) return true;
+
+    pmVisualInitWindow(&kapa, "ppmerge");
+    return true;
+}
+
+bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean) {
+
+    Graphdata graphdata;
+    float xline[2], yline[2];
+    
+    if (!pmVisualIsVisual()) return true;
+
+    if (!plotFlag) return true;
+
+    KapaInitGraph(&graphdata);
+
+    psVector *xAll = psVectorAlloc(pixels->n, PS_TYPE_F32);
+    psVector *xSub = psVectorAlloc(pixels->n, PS_TYPE_F32);
+    psVector *ySub = psVectorAlloc(pixels->n, PS_TYPE_F32);
+
+    // generate vectors of the unmasked values
+    int nSub = 0;
+    for (int j = 0; j < pixels->n; j++) {
+        xAll->data.F32[j] = j;
+        if (mask && mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) continue;
+        xSub->data.F32[nSub] = j;
+        ySub->data.F32[nSub] = pixels->data.F32[j];
+        nSub ++;
+    }
+    xSub->n = ySub->n = nSub;
+    xAll->n = pixels->n;
+        
+    xline[0] = 0;
+    xline[1] = pixels->n;
+    yline[0] = mean;
+    yline[1] = mean;
+
+    // plot the unmasked values
+    pmVisualScaleGraphdata (&graphdata, xAll, pixels, false);
+    KapaSetGraphData(kapa, &graphdata);
+    KapaSetLimits(kapa, &graphdata);
+    KapaClearPlots (kapa);
+
+    KapaSetFont (kapa, "courier", 14);
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "ordinate", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "pixel values", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName("black");
+    graphdata.style = 2;
+    graphdata.ptype = 2;
+    KapaPrepPlot  (kapa, xSub->n, &graphdata);
+    KapaPlotVector(kapa, xSub->n, xSub->data.F32, "x");
+    KapaPlotVector(kapa, xSub->n, ySub->data.F32, "y");
+
+    graphdata.color = KapaColorByName("red");
+    graphdata.style = 2;
+    graphdata.ptype = 7;
+    KapaPrepPlot  (kapa, xAll->n, &graphdata);
+    KapaPlotVector(kapa, xAll->n, xAll->data.F32, "x");
+    KapaPlotVector(kapa, xAll->n, pixels->data.F32, "y");
+
+    graphdata.color = KapaColorByName("blue");
+    graphdata.style = 0;
+    graphdata.ptype = 7;
+    KapaPrepPlot  (kapa, 2, &graphdata);
+    KapaPlotVector(kapa, 2, xline, "x");
+    KapaPlotVector(kapa, 2, yline, "y");
+
+    pmVisualAskUser (&plotFlag);
+    return true;
+}
+
+bool pmReadoutCombineVisualCleanup(void) {
+
+    if (!pmVisualIsVisual()) return true;
+    if (kapa == -1) return true;
+
+    KapaClose(kapa);
+    return true;
+}
+
+# else
+
+bool pmReadoutCombineVisualInit(void) { return true; }
+bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean) { return true; }
+bool pmReadoutCombineVisualCleanup(void) { return true; }
+
+# endif

branches/eam_branches/20090715/psModules/src/imcombine/pmReadoutCombine.h

@@ -47 +47 @@
                      );
 
+bool pmReadoutCombineVisualInit(void);
+bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean);
+bool pmReadoutCombineVisualCleanup(void);
+
 /// @}
 #endif

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionEquation.c

@@ -793 +793 @@
 
 #ifdef TESTING
-              // XXX double-check for NAN in data:
-                for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
-                    for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
-                        if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
-                            fprintf (stderr, "WARNING: NAN in matrix1\n");
-                        }
-                    }
-                }
-                for (int ix = 0; ix < stamp->vector1->n; ix++) {
-                    if (!isfinite(stamp->vector1->data.F64[ix])) {
-                        fprintf (stderr, "WARNING: NAN in vector1\n");
-                    }
-                }
+              // XXX double-check for NAN in data:
+                for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
+                    for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
+                        if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
+                            fprintf (stderr, "WARNING: NAN in matrix1\n");
+                        }
+                    }
+                }
+                for (int ix = 0; ix < stamp->vector1->n; ix++) {
+                    if (!isfinite(stamp->vector1->data.F64[ix])) {
+                        fprintf (stderr, "WARNING: NAN in vector1\n");
+                    }
+                }
 #endif
 
@@ -818 +818 @@
 
 #ifdef TESTING
-        for (int ix = 0; ix < sumVector->n; ix++) {
-            if (!isfinite(sumVector->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector1\n");
-            }
-        }
+        for (int ix = 0; ix < sumVector->n; ix++) {
+            if (!isfinite(sumVector->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
 #endif
 
@@ -828 +828 @@
 
 #ifdef TESTING
-        for (int ix = 0; ix < sumVector->n; ix++) {
-            if (!isfinite(sumVector->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector1\n");
-            }
-        }
+        for (int ix = 0; ix < sumVector->n; ix++) {
+            if (!isfinite(sumVector->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
         {
             psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
@@ -866 +866 @@
 
 #ifdef TESTING
-        // XXX double-check for NAN in data:
-        for (int ix = 0; ix < kernels->solution1->n; ix++) {
-            if (!isfinite(kernels->solution1->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector1\n");
-            }
-        }
+        // XXX double-check for NAN in data:
+        for (int ix = 0; ix < kernels->solution1->n; ix++) {
+            if (!isfinite(kernels->solution1->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
 #endif
 
@@ -1127 +1127 @@
                 source = stamp->image1;
                 convolutions = stamp->convolutions1;
+
+                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
+                // so that it is on the scale of image1.
+                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
+                                                          false); // Kernel image
+                if (!image) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
+                    return false;
+                }
+                double sumKernel = 0;   // Sum of kernel, for normalising residual
+                int size = kernels->size; // Half-size of kernel
+                int fullSize = 2 * size + 1; // Full size of kernel
+                for (int y = 0; y < fullSize; y++) {
+                    for (int x = 0; x < fullSize; x++) {
+                        sumKernel += image->data.F32[y][x];
+                    }
+                }
+                psFree(image);
+                devNorm = 1.0 / sumKernel / numPixels;
                 break;
               case PM_SUBTRACTION_MODE_2: