Changeset 24951 for branches/czw_branch/cleanup/psModules/src/imcombine

Timestamp:

Jul 30, 2009, 5:20:29 PM (17 years ago)

Author:

watersc1

Message:

Finished up my edits to the detrend cleanup, and some changes to my
copy of burntool and the pslib astrometry. Detrend cleanup has not
been tested yet. That's up next.

Location:

branches/czw_branch/cleanup

Files:

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

branches/czw_branch/cleanup/psModules/src/imcombine/pmReadoutCombine.c

@@ -173 +173 @@
     for (int i = 0; i < inputs->n; i++) {
         pmReadout *readout = inputs->data[i]; // Readout of interest
+        psAssert(readout, "readout was not defined");
+        if (!readout->process) continue;
         if (!readout->image) {
             psError(PS_ERR_UNEXPECTED_NULL, true, "Image data is missing for image %d.\n", i);
@@ -271 +273 @@
         for (int r = 0; r < inputs->n; r++) {
             pmReadout *readout = inputs->data[r]; // Input readout
+            if (!readout->process) continue; 
             psTrace("psModules.imcombine", 3, "Iterating input %d: %d --> %d, %d --> %d\n", r,
                     minInputCols - readout->col0, maxInputCols - readout->col0,
@@ -277 +280 @@
     }
     #endif
+
+    // set up windows for visualization (if selected)
+    pmReadoutCombineVisualInit();
 
     // Dereference output products
@@ -308 +314 @@
             for (int r = 0; r < inputs->n; r++) {
                 pmReadout *readout = inputs->data[r]; // Input readout
+                if (!readout->process) {
+                    maskData[r] = 1;
+                    continue;
+                }
                 int yIn = i - readout->row0; // y position on input readout
                 int xIn = j - readout->col0; // x position on input readout
@@ -390 +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;
@@ -427 +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/czw_branch/cleanup/psModules/src/imcombine/pmReadoutCombine.h

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

branches/czw_branch/cleanup/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: