Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/imcombine/pmSubtractionVisual.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionVisual.c (modified) (14 diffs)

branches/simtest_nebulous_branches/psModules/src/imcombine/pmSubtractionVisual.c

@@ -16 +16 @@
 
 #include "pmKapaPlots.h"
+#include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionEquation.h"
 
 #include "pmVisual.h"
@@ -32 +34 @@
 static bool plotLeastSquares     = true;
 static bool plotImage            = true;
+
 // variables to store plotting window indices
-static int kapa  = -1;
+static int kapa1 = -1;
 static int kapa2 = -1;
+static int kapa3 = -1;
 
 /** function prototypes*/
@@ -46 +50 @@
 bool pmSubtractionVisualClose(void)
 {
-    if(kapa != -1)
-        KiiClose(kapa);
-    if(kapa2 != -1)
-        KiiClose(kapa2);
+    if(kapa1 != -1) KiiClose(kapa1);
+    if(kapa2 != -1) KiiClose(kapa2);
     return true;
 }
@@ -60 +62 @@
 bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {
     if (!pmVisualIsVisual() || !plotConvKernels) return true;
-    if (!pmVisualInitWindow(&kapa, "ppSub:Images")) {
-        return false;
-    }
-    pmVisualScaleImage(kapa, convKernels, "Convolution_Kernels", 0, true);
+    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) {
+        return false;
+    }
+    pmVisualScaleImage(kapa1, convKernels, "Convolution_Kernels", 0, true);
     pmVisualAskUser(&plotConvKernels);
     return true;
@@ -74 +76 @@
 bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
     if (!pmVisualIsVisual() || !plotStamps) return true;
-    if (!pmVisualInitWindow (&kapa, "ppSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) {
         return false;
     }
@@ -134 +136 @@
         stampNum++;
     }
-    pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true);
+    pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);
 
     pmVisualAskUser(&plotStamps);
@@ -144 +146 @@
 
     if (!pmVisualIsVisual() || !plotLeastSquares) return true;
-    if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
         return false;
     }
@@ -157 +159 @@
         if (stamp == NULL) continue;
 
-        psImage *im = stamp->matrix1;
-        if (im == NULL) im = stamp->matrix2;
-        if (im == NULL) im = stamp->matrixX;
+        psImage *im = stamp->matrix;
         if (im == NULL) continue;
 
@@ -187 +187 @@
         if (stamp == NULL) continue;
 
-        psImage *im = stamp->matrix1;
-        if (im == NULL) im = stamp->matrix2;
-        if (im == NULL) im = stamp->matrixX;
+        psImage *im = stamp->matrix;
         if (im == NULL) continue;
 
@@ -197 +195 @@
 
     psImage *canvas32 = pmVisualImageToFloat(canvas);
-    pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true);
+    pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
 
     pmVisualAskUser(&plotLeastSquares);
@@ -207 +205 @@
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
     if (!pmVisualIsVisual() || !plotImage) return true;
-    if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
-        return false;
-    }
-
-    pmVisualScaleImage(kapa, image, "Image", 0, true);
-    pmVisualScaleImage(kapa, ref, "Reference", 1, true);
-    pmVisualScaleImage(kapa, sub, "Subtraction", 2, true);
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
+        return false;
+    }
+
+    pmVisualScaleImage(kapa1, image, "Image", 0, true);
+    pmVisualScaleImage(kapa1, ref, "Reference", 1, true);
+    pmVisualScaleImage(kapa1, sub, "Subtraction", 2, true);
+    pmVisualAskUser(&plotImage);
+    return true;
+}
+
+bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) {
+
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
+        return false;
+    }
+
+    // get the kernel sizes
+    int footprint = kernels->size;
+
+    // output image is a grid of NXsub by NYsub sub-images
+    int NXsub = sqrt(kernels->num);
+    int NYsub = kernels->num / NXsub;
+    if (kernels->num % NXsub) NYsub++;
+
+    int NXpix = NXsub * (2*footprint + 1 + 3);
+    int NYpix = NYsub * (2*footprint + 1 + 3);
+
+    psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
+    psImageInit (output, 0.0);
+
+    for (int i = 0; i < kernels->num; i++) {
+        pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i];
+        psKernel *kernel = preCalc->kernel;
+
+        int xSub = i % NXsub;
+        int ySub = i / NXsub;
+
+        int xPix = xSub * (2*footprint + 1 + 3) + footprint;
+        int yPix = ySub * (2*footprint + 1 + 3) + footprint;
+
+        double sum = 0.0;
+        for (int y = -footprint; y <= footprint; y++) {
+            for (int x = -footprint; x <= footprint; x++) {
+                output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                sum += kernel->kernel[y][x];
+            }
+        }
+        fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+    }                                                    
+        
+    pmVisualScaleImage(kapa1, output, "Image", 0, true);
+    pmVisualAskUser(&plotImage);
+    return true;
+}
+
+bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps) {
+
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) {
+        return false;
+    }
+
+    // get the kernel sizes
+    int footprint = stamps->footprint;
+
+    // choose the brightest stamp
+    pmSubtractionStamp *maxStamp = NULL;
+    float maxFlux = NAN;
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i];
+        if (!isfinite(stamp->flux)) continue;
+        if (!stamp->convolutions1 && !stamp->convolutions2) continue;
+        if (!maxStamp) {
+            maxFlux = stamp->flux;
+            maxStamp = stamp;
+            continue;
+        }
+        if (stamp->flux > maxFlux) {
+            maxFlux = stamp->flux;
+            maxStamp = stamp;
+        }
+    }
+
+    if (!isfinite(maxStamp->flux)) {
+        fprintf (stderr, "no valid stamps?\n");
+    }
+
+    int nKernels = 0;
+
+    if (maxStamp->convolutions1) {
+        // output image is a grid of NXsub by NYsub sub-images
+        nKernels = maxStamp->convolutions1->n;
+        int NXsub = sqrt(nKernels);
+        int NYsub = nKernels / NXsub;
+        if (nKernels % NXsub) NYsub++;
+
+        int NXpix = NXsub * (2*footprint + 1 + 3);
+        int NYpix = NYsub * (2*footprint + 1 + 3);
+
+        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
+        psImageInit (output, 0.0);
+
+        for (int i = 0; i < nKernels; i++) {
+            psKernel *kernel = maxStamp->convolutions1->data[i];
+            
+            int xSub = i % NXsub;
+            int ySub = i / NXsub;
+            
+            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
+            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
+            
+            double sum = 0.0;
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                    sum += kernel->kernel[y][x];
+                }
+            }
+            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+        }               
+        pmVisualScaleImage(kapa2, output, "Image", 0, true);
+    }                                    
+        
+    if (maxStamp->convolutions2) {
+        // output image is a grid of NXsub by NYsub sub-images
+        nKernels = maxStamp->convolutions2->n;
+        int NXsub = sqrt(nKernels);
+        int NYsub = nKernels / NXsub;
+        if (nKernels % NXsub) NYsub++;
+
+        int NXpix = NXsub * (2*footprint + 1 + 3);
+        int NYpix = NYsub * (2*footprint + 1 + 3);
+
+        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
+        psImageInit (output, 0.0);
+
+        for (int i = 0; i < nKernels; i++) {
+            psKernel *kernel = maxStamp->convolutions2->data[i];
+            
+            int xSub = i % NXsub;
+            int ySub = i / NXsub;
+            
+            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
+            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
+            
+            double sum = 0.0;
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                    sum += kernel->kernel[y][x];
+                }
+            }
+            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+        }               
+        pmVisualScaleImage(kapa2, output, "Image", 1, true);
+    }                                    
+        
     pmVisualAskUser(&plotImage);
     return true;
@@ -240 +390 @@
 
         overlay[Noverlay].type = KII_OVERLAY_BOX;
+        if ((stamp->x < 1.0) && (stamp->y < 1.0)) {
+            // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
+            continue;
+        }
+        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
+            // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
+            continue;
+        }
         overlay[Noverlay].x = stamp->x;
         overlay[Noverlay].y = stamp->y;
@@ -255 +413 @@
 }
 
+static int footprint = 0;
+static int NX = 0;
+static int NY = 0;
+static psImage *sourceImage      = NULL;
+static psImage *targetImage      = NULL;
+static psImage *residualImage    = NULL;
+static psImage *fresidualImage   = NULL;
+static psImage *differenceImage  = NULL;
+static psImage *convolutionImage = NULL;
+
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
+
+    if (!pmVisualIsVisual()) return true;
+
+    // generate 4 storage images large enough to hold the N stamps:
+
+    footprint = stamps->footprint;
+
+    float NXf = sqrt(stamps->num);
+    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
+    
+    float NYf = stamps->num / NX;
+    NY = (int) NYf == NY ? NYf : NYf + 1.0;
+
+    int NXpix = (2*footprint + 1) * NX;
+    NXpix += (NX > 1) ? 3 * NX : 0;
+
+    int NYpix = (2*footprint + 1) * NY;
+    NYpix += (NY > 1) ? 3 * NY : 0;
+
+    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    
+    psImageInit (sourceImage,      0.0);
+    psImageInit (targetImage,      0.0);
+    psImageInit (residualImage,    0.0);
+    psImageInit (fresidualImage,   0.0);
+    psImageInit (differenceImage,  0.0);
+    psImageInit (convolutionImage, 0.0);
+
+    return true;
+}
+
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
+
+    if (!pmVisualIsVisual()) return true;
+
+    double sum;
+
+    int NXoff = index % NX;
+    int NYoff = index / NX;
+
+    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
+    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
+
+    // insert the (target) kernel into the (target) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
+            sum += targetImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (source) kernel into the (source) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
+            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (convolution) kernel into the (convolution) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x];
+            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    
+    // insert the (difference) kernel into the (difference) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
+            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (residual) kernel into the (residual) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
+            sum += residualImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (fresidual) kernel into the (fresidual) image:
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
+        }
+    }
+    return true;
+}
+
+bool pmSubtractionVisualShowFit(double norm) {
+
+    // for testing, dump the residual image and exit
+    if (0) {
+        psMetadata *header = psMetadataAlloc();
+        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
+        psFits *fits = psFitsOpen("resid.fits", "w");
+        psFitsWriteImage(fits, header, residualImage, 0, NULL);
+        psFitsClose(fits);
+        // exit (0);
+    }
+
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
+    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
+
+    KiiEraseOverlay (kapa1, "red");
+    KiiEraseOverlay (kapa2, "red");
+
+    pmVisualScaleImage(kapa1, targetImage, "Target Stamps", 0, true);
+    pmVisualScaleImage(kapa1, sourceImage, "Source Stamps", 1, true);
+    pmVisualScaleImage(kapa1, convolutionImage, "Convolution Stamps", 2, true);
+    KiiCenter (kapa1, 0.5*targetImage->numCols, 0.5*targetImage->numRows, 1);
+
+    pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true);
+    pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);
+
+    if (1) {
+        KiiImage image;
+        KapaImageData data;
+        Coords coords;
+        strcpy (coords.ctype, "RA---TAN");
+
+        image.data2d = residualImage->data.F32;
+        image.Nx = residualImage->numCols;
+        image.Ny = residualImage->numRows;
+        strcpy (data.name, "Residual Stamps");
+        strcpy (data.file, "Residual Stamps");
+
+        data.zero  = -300.0;
+        data.range = +600.0;
+        data.logflux = 0;
+
+        KiiSetChannel (kapa2, 1);
+        KiiNewPicture2D (kapa2, &image, &data, &coords);
+    } else {
+        pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true);
+    }
+
+    KiiCenter (kapa2, 0.5*residualImage->numCols, 0.5*residualImage->numRows, 1);
+
+    pmVisualAskUser(NULL);
+
+    psFree(targetImage);
+    psFree(sourceImage);
+    psFree(convolutionImage);
+    psFree(differenceImage);
+    psFree(residualImage);
+    psFree(fresidualImage);
+
+    targetImage = NULL;
+    sourceImage = NULL;
+    convolutionImage = NULL;
+    differenceImage = NULL;
+    residualImage = NULL;
+    fresidualImage = NULL;
+
+    return true;
+}
+
+bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels) {
+
+    Graphdata graphdata;
+
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;
+
+    KapaClearSections (kapa3);
+    KapaInitGraph (&graphdata);
+
+    psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
+    psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
+
+    graphdata.xmin = -1.0;
+    graphdata.xmax = kernels->num + 1.0;
+    graphdata.ymin = +32.0;
+    graphdata.ymax = -32.0;
+
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0);
+
+    // construct the plot vectors
+    for (int i = 0; i < kernels->num; i++) {
+        x->data.F32[i] = i;
+        y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
+        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
+        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
+    }
+    x->n = y->n = kernels->num;
+
+    float range;
+    range = graphdata.xmax - graphdata.xmin;
+    graphdata.xmax += 0.05*range;
+    graphdata.xmin -= 0.05*range;
+    range = graphdata.ymax - graphdata.ymin;
+    graphdata.ymax += 0.05*range;
+    graphdata.ymin -= 0.05*range;
+
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "kernel number", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "coeff", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    KapaPrepPlot   (kapa3, x->n, &graphdata);
+    KapaPlotVector (kapa3, x->n, x->data.F32, "x");
+    KapaPlotVector (kapa3, x->n, y->data.F32, "y");
+
+    psFree (x);
+    psFree (y);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
 #else
 bool pmSubtractionVisualClose(void) {return true;}
@@ -261 +667 @@
 bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps) {return true;}
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {return true;}
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {return true;}
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {return true;}
+bool pmSubtractionVisualShowFit() {return true;}
+bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
 #endif