IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 13, 2011, 12:19:53 PM (15 years ago)
Author:
eugene
Message:

some code reorganzation (move all types into pmSubtractionTypes.h); add stage to check on different mode and order options, choosing best based on chisq of subtraction; careful with the window and the kernel sizes: if the measured kron radius is too large, allow the window to grow as needed; use the 1st radial moment measured on the stamps to set the kernel scaling; remove some cruft from the code; calling program now passes in kernel scaling options to be used when 1st radial moment is measured; new function to update the kernel description string used for kernel I/O; update the kernel description after the spatial order and direction is chosen; need to carry the kernel fwhms and spatial order to allow for update; calculate the psf-match solution errors; psf solution now uses the weights to get sensible chisq values, window is used to calculate the moments; penalty scale is now adjusted to be consistent with sensible reduced chisq; improved visualizations; use range-limiter in SVD of 1e-10; calculate chisq and moments for the solution and use in the evaluation; split out the stamp selection / convolution steps; reject sources after fitting a model to the chisq assuming non-zero systematic error

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtractionVisual.c

    r29600 r30622  
    1616
    1717#include "pmKapaPlots.h"
     18#include "pmFPA.h"
     19#include "pmSubtractionTypes.h"
    1820#include "pmSubtraction.h"
    1921#include "pmSubtractionStamps.h"
    2022#include "pmSubtractionEquation.h"
    2123#include "pmSubtractionKernels.h"
     24#include "pmSubtractionVisual.h"
    2225
    2326#include "pmVisual.h"
    … …  
    210213        psImageOverlaySection(canvas, im, x0, y0, "=");
    211214        stampNum++;
     215
     216        // renormalize the section
     217        float maxValue = 0;
     218        for (int iy = 0; iy < im->numRows; iy++) {
     219            for (int ix = 0; ix < im->numCols; ix++) {
     220                maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]);
     221            }
     222        }
     223        if (maxValue == 0.0) continue;
     224        for (int iy = 0; iy < im->numRows; iy++) {
     225            for (int ix = 0; ix < im->numCols; ix++) {
     226                canvas->data.F64[y0 + iy][x0 + ix] /= maxValue;
     227            }
     228        }
    212229    }
    213230
    … …  
    231248}
    232249
     250/** Plot the least-squares matrix of each stamp */
     251bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed) {
     252
     253    if (!pmVisualTestLevel("ppsub.chisq", 1)) return true;
     254
     255    if (!plotLeastSquares) return true;
     256
     257    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
     258        return false;
     259    }
     260
     261    psImage *matrixNorm = psImageCopy(NULL, matrixIn, PS_TYPE_F64);
     262    {
     263        // renormalize the matrix
     264        float maxValue = 0;
     265        for (int iy = 0; iy < matrixNorm->numRows; iy++) {
     266            for (int ix = 0; ix < matrixNorm->numCols; ix++) {
     267                maxValue = PS_MAX(maxValue, matrixNorm->data.F64[iy][ix]);
     268            }
     269        }
     270        for (int iy = 0; iy < matrixNorm->numRows; iy++) {
     271            for (int ix = 0; ix < matrixNorm->numCols; ix++) {
     272                matrixNorm->data.F64[iy][ix] /= maxValue;
     273            }
     274        }
     275    }
     276
     277    // Find the stamp size
     278    int imageMax = -1;
     279    int numStamps = 0;
     280    psElemType type = PS_TYPE_F64;
     281
     282    for(int i = 0; i < stamps->num; i++) {
     283        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     284        if (stamp == NULL) continue;
     285
     286        psImage *im = stamp->matrix;
     287        if (im == NULL) continue;
     288
     289        imageMax = PS_MAX(imageMax, im->numCols);
     290        imageMax = PS_MAX(imageMax, im->numRows);
     291        numStamps++;
     292        type = im->type.type;
     293    }
     294    if (imageMax == -1) return false;
     295
     296    int border = 15;
     297    imageMax += border;
     298    int tileRowCount = (int) ceil(sqrt(numStamps));
     299    int canvasX = tileRowCount * (imageMax);
     300    int canvasY = tileRowCount * (imageMax);
     301    psImage *canvas = psImageAlloc (canvasX, canvasY, type);
     302    psImageInit (canvas, NAN);
     303
     304    // overlay the images
     305    int stampNum = 0;
     306    int stampListNum = 0;
     307    while (stampNum < numStamps) {
     308        int x0 = (imageMax) * (stampNum % tileRowCount);
     309        int y0 = (imageMax) * (stampNum / tileRowCount);
     310
     311        pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++];
     312        if (stamp == NULL) continue;
     313
     314        psImage *im = stamp->matrix;
     315        if (im == NULL) continue;
     316
     317        stampNum++;
     318
     319        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     320
     321        // renormalize the section
     322        float maxValue = 0;
     323        for (int iy = 0; iy < im->numRows; iy++) {
     324            for (int ix = 0; ix < im->numCols; ix++) {
     325                maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]);
     326            }
     327        }
     328        if (maxValue == 0.0) continue;
     329        for (int iy = 0; iy < im->numRows; iy++) {
     330            for (int ix = 0; ix < im->numCols; ix++) {
     331                canvas->data.F64[y0 + iy][x0 + ix] = (im->data.F64[iy][ix] / maxValue - matrixNorm->data.F64[iy][ix]) / matrixNorm->data.F64[iy][ix];
     332            }
     333        }
     334    }
     335
     336    psImage *canvas32 = pmVisualImageToFloat(canvas);
     337    pmVisualRangeImage(kapa2, canvas32, "Least_Squares", 0, -100.0, 100.0);
     338
     339    if (0) {
     340        static int count = 0;
     341        char filename[64];
     342        sprintf (filename, "chisq.%02d.fits", count);
     343        count ++;
     344        psFits *fits = psFitsOpen (filename, "w");
     345        psFitsWriteImage (fits, NULL, canvas32, 0, NULL);
     346        psFitsClose (fits);
     347    }
     348
     349    pmVisualAskUser(&plotLeastSquares);
     350    psFree(canvas);
     351    psFree(canvas32);
     352    return true;
     353}
     354
    233355bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
    234356
    … …  
    304426    }
    305427
    306     // XXX clear the overlay(s) (red at least!)
     428    // clear the overlay (red at least!)
     429    KiiEraseOverlay (kapa2, "red");
    307430
    308431    // get the kernel sizes
    … …  
    337460    int nKernels = 0;
    338461
     462    // paste in the kernel images, scaled by sum2
    339463    if (maxStamp->convolutions1) {
    340464        // output image is a grid of NXsub by NYsub sub-images
    … …  
    359483            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
    360484           
    361             double sum = 0.0;
    362485            double sum2 = 0.0;
    363486            for (int y = -footprint; y <= footprint; y++) {
    364487                for (int x = -footprint; x <= footprint; x++) {
    365                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
    366                     sum += kernel->kernel[y][x];
    367488                    sum2 += PS_SQR(kernel->kernel[y][x]);
    368489                }
    369490            }
    370             // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
     491            float scale = sqrt(sum2) / PS_SQR(2*footprint + 1);
     492            for (int y = -footprint; y <= footprint; y++) {
     493                for (int x = -footprint; x <= footprint; x++) {
     494                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale;
     495                }
     496            }
    371497        }               
    372498        pmVisualScaleImage(kapa2, output, "Image", 0, true);
    … …  
    401527            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
    402528           
    403             double sum = 0.0;
    404529            double sum2 = 0.0;
    405530            for (int y = -footprint; y <= footprint; y++) {
    406531                for (int x = -footprint; x <= footprint; x++) {
    407                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
    408                     sum += kernel->kernel[y][x];
    409532                    sum2 += PS_SQR(kernel->kernel[y][x]);
    410533                }
    411534            }
    412             // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
     535            float scale = sqrt(sum2) / PS_SQR(2*footprint + 1);
     536            for (int y = -footprint; y <= footprint; y++) {
     537                for (int x = -footprint; x <= footprint; x++) {
     538                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale;
     539                }
     540            }
    413541        }               
    414542        pmVisualScaleImage(kapa2, output, "Image", 1, true);
    … …  
    468596    KiiLoadOverlay (kapa2, overlay, Noverlay, "red");
    469597    FREE (overlay);
    470     return true;
    471 }
    472 
    473 static int footprint = 0;
    474 static int NX = 0;
    475 static int NY = 0;
    476 static psImage *sourceImage      = NULL;
    477 static psImage *targetImage      = NULL;
    478 static psImage *residualImage    = NULL;
    479 static psImage *fresidualImage   = NULL;
    480 static psImage *differenceImage  = NULL;
    481 static psImage *convolutionImage = NULL;
    482 
    483 bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
    484 
    485     if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
    486 
    487     // generate 4 storage images large enough to hold the N stamps:
    488 
    489     footprint = stamps->footprint;
    490 
    491     float NXf = sqrt(stamps->num);
    492     NX = (int) NXf == NXf ? NXf : NXf + 1.0;
    493    
    494     float NYf = stamps->num / NX;
    495     NY = (int) NYf == NY ? NYf : NYf + 1.0;
    496 
    497     int NXpix = (2*footprint + 1) * NX;
    498     NXpix += (NX > 1) ? 3 * NX : 0;
    499 
    500     int NYpix = (2*footprint + 1) * NY;
    501     NYpix += (NY > 1) ? 3 * NY : 0;
    502 
    503     sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    504     targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    505     residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    506     fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    507     differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    508     convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    509    
    510     psImageInit (sourceImage,      0.0);
    511     psImageInit (targetImage,      0.0);
    512     psImageInit (residualImage,    0.0);
    513     psImageInit (fresidualImage,   0.0);
    514     psImageInit (differenceImage,  0.0);
    515     psImageInit (convolutionImage, 0.0);
    516 
    517     return true;
    518 }
    519 
    520 bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
    521 
    522     if (!pmVisualTestLevel("ppsub.stamp", 1)) return true;
    523 
    524     double sum;
    525 
    526     int NXoff = index % NX;
    527     int NYoff = index / NX;
    528 
    529     int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
    530     int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
    531 
    532     // insert the (target) kernel into the (target) image:
    533     sum = 0.0;
    534     for (int y = -footprint; y <= footprint; y++) {
    535         for (int x = -footprint; x <= footprint; x++) {
    536             targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
    537             sum += targetImage->data.F32[y + NYpix][x + NXpix];
    538         }
    539     }
    540     targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    541 
    542     // insert the (source) kernel into the (source) image:
    543     sum = 0.0;
    544     for (int y = -footprint; y <= footprint; y++) {
    545         for (int x = -footprint; x <= footprint; x++) {
    546             sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
    547             sum += sourceImage->data.F32[y + NYpix][x + NXpix];
    548         }
    549     }
    550     sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    551 
    552     // insert the (convolution) kernel into the (convolution) image:
    553     sum = 0.0;
    554     for (int y = -footprint; y <= footprint; y++) {
    555         for (int x = -footprint; x <= footprint; x++) {
    556             convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x];
    557             sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
    558         }
    559     }
    560     convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    561    
    562     // insert the (difference) kernel into the (difference) image:
    563     sum = 0.0;
    564     for (int y = -footprint; y <= footprint; y++) {
    565         for (int x = -footprint; x <= footprint; x++) {
    566             differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
    567             sum += differenceImage->data.F32[y + NYpix][x + NXpix];
    568         }
    569     }
    570     differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    571 
    572     // insert the (residual) kernel into the (residual) image:
    573     sum = 0.0;
    574     for (int y = -footprint; y <= footprint; y++) {
    575         for (int x = -footprint; x <= footprint; x++) {
    576             residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
    577             sum += residualImage->data.F32[y + NYpix][x + NXpix];
    578         }
    579     }
    580     residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    581 
    582     // insert the (fresidual) kernel into the (fresidual) image:
    583     for (int y = -footprint; y <= footprint; y++) {
    584         for (int x = -footprint; x <= footprint; x++) {
    585             fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
    586         }
    587     }
    588598    return true;
    589599}
    … …  
    611621}
    612622
    613 bool pmSubtractionVisualShowFit(double norm) {
    614 
    615     // for testing, dump the residual image and exit
    616     if (0) {
    617         psMetadata *header = psMetadataAlloc();
    618         psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
    619         psFits *fits = psFitsOpen("resid.fits", "w");
    620         psFitsWriteImage(fits, header, residualImage, 0, NULL);
    621         psFitsClose(fits);
    622         // exit (0);
    623     }
     623static int footprint = 0;
     624static int NX = 0;
     625static int NY = 0;
     626static psImage *sourceImage      = NULL;
     627static psImage *targetImage      = NULL;
     628static psImage *residualImage    = NULL;
     629static psImage *fresidualImage   = NULL;
     630static psImage *differenceImage  = NULL;
     631static psImage *convolutionImage = NULL;
     632
     633bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) {
    624634
    625635    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
    … …  
    627637    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
    628638    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
     639
     640    // set up holding images for the visualization
     641    pmSubtractionVisualShowFitInit (stamps);
     642
     643    int numKernels = kernels->num;      // Number of kernels
     644
     645    psImage *polyValues = NULL;         // Polynomial values
     646    psKernel *residual = psKernelAlloc(-stamps->footprint, stamps->footprint, -stamps->footprint, stamps->footprint); // Residual image
     647
     648    double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     649
     650    for (int i = 0; i < stamps->num; i++) {
     651        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
     652        if (stamp->status != PM_SUBTRACTION_STAMP_USED) { continue; }
     653
     654        // Calculate coefficients of the kernel basis functions
     655        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     656        double background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Difference in background
     657
     658        psImageInit(residual->image, 0.0);
     659
     660        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
     661            psKernel *target;           // Target postage stamp
     662            psKernel *source;           // Source postage stamp
     663            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
     664            switch (kernels->mode) {
     665              case PM_SUBTRACTION_MODE_1:
     666                target = stamp->image2;
     667                source = stamp->image1;
     668                convolutions = stamp->convolutions1;
     669                break;
     670              case PM_SUBTRACTION_MODE_2:
     671                target = stamp->image1;
     672                source = stamp->image2;
     673                convolutions = stamp->convolutions2;
     674                break;
     675              default:
     676                psAbort("Unsupported subtraction mode: %x", kernels->mode);
     677            }
     678
     679            for (int j = 0; j < numKernels; j++) {
     680                psKernel *convolution = convolutions->data[j]; // Convolution
     681                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
     682                for (int y = - footprint; y <= footprint; y++) {
     683                    for (int x = - footprint; x <= footprint; x++) {
     684                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
     685                    }
     686                }
     687            }
     688            // visualize the target, source, convolution and residual
     689            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
     690        } else {
     691            // Dual convolution
     692            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
     693            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
     694            psKernel *image1 = stamp->image1; // The first image
     695            psKernel *image2 = stamp->image2; // The second image
     696
     697            for (int j = 0; j < numKernels; j++) {
     698                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
     699                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
     700                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
     701                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
     702
     703                for (int y = - footprint; y <= footprint; y++) {
     704                    for (int x = - footprint; x <= footprint; x++) {
     705                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
     706                    }
     707                }
     708            }
     709            // visualize the target, source, convolution and residual
     710            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
     711        }
     712    }
     713    pmSubtractionVisualShowFitImage(norm);
     714
     715    return true;
     716}
     717
     718// generate 4 storage images large enough to hold the stamps:
     719bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
     720
     721    footprint = stamps->footprint;
     722
     723    float NXf = sqrt(stamps->num);
     724    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
     725   
     726    float NYf = stamps->num / NX;
     727    NY = (int) NYf == NY ? NYf : NYf + 1.0;
     728
     729    int NXpix = (2*footprint + 1) * NX;
     730    NXpix += (NX > 1) ? 3 * NX : 0;
     731
     732    int NYpix = (2*footprint + 1) * NY;
     733    NYpix += (NY > 1) ? 3 * NY : 0;
     734
     735    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     736    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     737    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     738    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     739    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     740    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     741   
     742    psImageInit (sourceImage,      0.0);
     743    psImageInit (targetImage,      0.0);
     744    psImageInit (residualImage,    0.0);
     745    psImageInit (fresidualImage,   0.0);
     746    psImageInit (differenceImage,  0.0);
     747    psImageInit (convolutionImage, 0.0);
     748
     749    return true;
     750}
     751
     752bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
     753
     754    double sum;
     755
     756    int NXoff = index % NX;
     757    int NYoff = index / NX;
     758
     759    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
     760    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
     761
     762    // insert the (target) kernel into the (target) image:
     763    sum = 0.0;
     764    for (int y = -footprint; y <= footprint; y++) {
     765        for (int x = -footprint; x <= footprint; x++) {
     766            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
     767            sum += targetImage->data.F32[y + NYpix][x + NXpix];
     768        }
     769    }
     770    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     771
     772    // insert the (source) kernel into the (source) image:
     773    sum = 0.0;
     774    for (int y = -footprint; y <= footprint; y++) {
     775        for (int x = -footprint; x <= footprint; x++) {
     776            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
     777            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
     778        }
     779    }
     780    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     781
     782    // insert the (convolution) kernel into the (convolution) image:
     783    sum = 0.0;
     784    for (int y = -footprint; y <= footprint; y++) {
     785        for (int x = -footprint; x <= footprint; x++) {
     786            convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x];
     787            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
     788        }
     789    }
     790    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     791   
     792    // insert the (difference) kernel into the (difference) image:
     793    sum = 0.0;
     794    for (int y = -footprint; y <= footprint; y++) {
     795        for (int x = -footprint; x <= footprint; x++) {
     796            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
     797            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
     798        }
     799    }
     800    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     801
     802    // insert the (residual) kernel into the (residual) image:
     803    sum = 0.0;
     804    for (int y = -footprint; y <= footprint; y++) {
     805        for (int x = -footprint; x <= footprint; x++) {
     806            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
     807            sum += residualImage->data.F32[y + NYpix][x + NXpix];
     808        }
     809    }
     810    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     811
     812    // insert the (fresidual) kernel into the (fresidual) image:
     813    for (int y = -footprint; y <= footprint; y++) {
     814        for (int x = -footprint; x <= footprint; x++) {
     815            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
     816        }
     817    }
     818    return true;
     819}
     820
     821bool pmSubtractionVisualShowFitImage(double norm) {
    629822
    630823    KiiEraseOverlay (kapa1, "red");
    … …  
    673866    psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
    674867    psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
     868    psVector *dy = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
    675869
    676870    graphdata.xmin = -1.0;
    … …  
    685879        x->data.F32[i] = i;
    686880        y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
     881        dy->data.F32[i] = kernels->solution1err->data.F64[i];
    687882        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
    688883        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
    689884    }
    690     x->n = y->n = kernels->num;
     885    x->n = y->n = dy->n = kernels->num;
    691886
    692887    float range;
    … …  
    709904    graphdata.size = 0.5;
    710905    graphdata.style = 2;
     906    graphdata.etype |= 0x01;
    711907
    712908    KapaPrepPlot   (kapa3, x->n, &graphdata);
    713909    KapaPlotVector (kapa3, x->n, x->data.F32, "x");
    714910    KapaPlotVector (kapa3, x->n, y->data.F32, "y");
     911    KapaPlotVector (kapa3, x->n, dy->data.F32, "dym");
     912    KapaPlotVector (kapa3, x->n, dy->data.F32, "dyp");
    715913
    716914    psFree (x);
    717915    psFree (y);
     916    psFree (dy);
    718917    psFree (polyValues);
     918
     919    pmVisualAskUser(NULL);
     920    return true;
     921}
     922
     923// plot log(flux) vs log(chisq), log(flux) vs log(moments), log(chisq) vs log(moments)
     924bool pmSubtractionVisualPlotChisqAndMoments(psVector *fluxes, psVector *chisq, psVector *moments) {
     925
     926    Graphdata graphdata;
     927    KapaSection section;
     928
     929    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
     930
     931    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;
     932
     933    KapaClearSections (kapa3);
     934    KapaInitGraph (&graphdata);
     935    KiiResize(kapa3, 1500, 500);
     936
     937    psVector *lchi = psVectorAlloc (fluxes->n, PS_TYPE_F32);
     938    psVector *lflx = psVectorAlloc (fluxes->n, PS_TYPE_F32);
     939    psVector *lMxx = psVectorAlloc (fluxes->n, PS_TYPE_F32);
     940
     941    // construct the plot vectors
     942    for (int i = 0; i < fluxes->n; i++) {
     943        lchi->data.F32[i] = log10(chisq->data.F32[i]);
     944        lflx->data.F32[i] = log10(fluxes->data.F32[i]);
     945        lMxx->data.F32[i] = log10(moments->data.F32[i]);
     946    }
     947
     948    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
     949
     950    // section 1: lflux vs lchi
     951    section.dx = 0.33;
     952    section.dy = 1.00;
     953    section.x  = 0.00;
     954    section.y  = 0.00;
     955    section.name = psStringCopy ("flux.v.chi");
     956    KapaSetSection (kapa3, &section);
     957    psFree (section.name);
     958
     959    graphdata.color = KapaColorByName ("black");
     960    pmVisualScaleGraphdata(&graphdata, lflx, lchi, false);
     961    KapaSetLimits (kapa3, &graphdata);
     962
     963    KapaSetFont (kapa3, "helvetica", 14);
     964    KapaBox (kapa3, &graphdata);
     965    KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM);
     966    KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_YM);
     967
     968    graphdata.color = KapaColorByName ("black");
     969    graphdata.ptype = 2;
     970    graphdata.size = 0.5;
     971    graphdata.style = 2;
     972
     973    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
     974    KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x");
     975    KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "y");
     976
     977    // section 2: lflux vs lMxx
     978    section.dx = 0.33;
     979    section.dy = 1.00;
     980    section.x  = 0.33;
     981    section.y  = 0.00;
     982    section.name = psStringCopy ("flux.v.mom");
     983    KapaSetSection (kapa3, &section);
     984    psFree (section.name);
     985
     986    graphdata.color = KapaColorByName ("black");
     987    pmVisualScaleGraphdata(&graphdata, lflx, lMxx, false);
     988    KapaSetLimits (kapa3, &graphdata);
     989
     990    KapaSetFont (kapa3, "helvetica", 14);
     991    KapaBox (kapa3, &graphdata);
     992    KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM);
     993    KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM);
     994
     995    graphdata.color = KapaColorByName ("black");
     996    graphdata.ptype = 2;
     997    graphdata.size = 0.5;
     998    graphdata.style = 2;
     999
     1000    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
     1001    KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x");
     1002    KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y");
     1003
     1004    // section 1: lflux vs lchi
     1005    section.dx = 0.33;
     1006    section.dy = 1.00;
     1007    section.x  = 0.66;
     1008    section.y  = 0.00;
     1009    section.name = psStringCopy ("chi.v.mom");
     1010    KapaSetSection (kapa3, &section);
     1011    psFree (section.name);
     1012
     1013    graphdata.color = KapaColorByName ("black");
     1014    pmVisualScaleGraphdata(&graphdata, lchi, lMxx, false);
     1015    KapaSetLimits (kapa3, &graphdata);
     1016
     1017    KapaSetFont (kapa3, "helvetica", 14);
     1018    KapaBox (kapa3, &graphdata);
     1019    KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_XM);
     1020    KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM);
     1021
     1022    graphdata.color = KapaColorByName ("black");
     1023    graphdata.ptype = 2;
     1024    graphdata.size = 0.5;
     1025    graphdata.style = 2;
     1026
     1027    KapaPrepPlot   (kapa3, lflx->n, &graphdata);
     1028    KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "x");
     1029    KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y");
    7191030
    7201031    pmVisualAskUser(NULL);
Note: See TracChangeset for help on using the changeset viewer.