Changeset 26573

Timestamp:

Jan 12, 2010, 3:53:09 PM (16 years ago)

Author:

eugene

Message:

sign has changed on convolution coefficients; require -visual as normal

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (modified) (16 diffs)

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c

@@ -238 +238 @@
 
     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);
-    }
-
+        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);
@@ -276 +276 @@
     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;
-        }
+        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");
+        fprintf (stderr, "no valid stamps?\n");
     }
 
@@ -297 +297 @@
 
     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++) {
+        // 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);
-    }
-
+            
+            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++) {
+        // 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);
-    }
-
+            
+            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;
@@ -390 +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;
-        }
+        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;
@@ -425 +425 @@
 bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
 
-    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualIsVisual()) return true;
 
     // generate 4 storage images large enough to hold the N stamps:
@@ -433 +433 @@
     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;
@@ -449 +449 @@
     differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-
+    
     psImageInit (sourceImage,      0.0);
     psImageInit (targetImage,      0.0);
@@ -462 +462 @@
 bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
 
-    // if (!pmVisualIsVisual()) return true;
+    if (!pmVisualIsVisual()) return true;
 
     double sum;
@@ -475 +475 @@
     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];
-        }
+        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;
@@ -485 +485 @@
     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];
-        }
+        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;
@@ -495 +495 @@
     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];
-        }
+        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];
-        }
+        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;
@@ -515 +515 @@
     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];
-        }
+        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;
@@ -524 +524 @@
     // 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));
-        }
+        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;
@@ -533 +533 @@
 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;
-
-    // XXX a dumb test : dump the residual image and exit
-    {
-        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 (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
     if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
@@ -558 +557 @@
     pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true);
     pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);
-    pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, 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);
 
@@ -603 +624 @@
     for (int i = 0; i < kernels->num; i++) {
         x->data.F32[i] = i;
-        y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
+        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]);