IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 10, 2010, 7:34:39 PM (16 years ago)
Author:
eugene
Message:

updates from eam_branches/20091201 (substantially changes to the kernel matching code; updates to pmDetection container, added pmPattern, etc)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r23487 r26893  
    1616
    1717#include "pmKapaPlots.h"
     18#include "pmSubtraction.h"
    1819#include "pmSubtractionStamps.h"
     20#include "pmSubtractionEquation.h"
    1921
    2022#include "pmVisual.h"
    … …  
    3234static bool plotLeastSquares     = true;
    3335static bool plotImage            = true;
     36
    3437// variables to store plotting window indices
    35 static int kapa  = -1;
     38static int kapa1 = -1;
    3639static int kapa2 = -1;
     40static int kapa3 = -1;
    3741
    3842/** function prototypes*/
    … …  
    4650bool pmSubtractionVisualClose(void)
    4751{
    48     if(kapa != -1)
    49         KiiClose(kapa);
    50     if(kapa2 != -1)
    51         KiiClose(kapa2);
     52    if(kapa1 != -1) KiiClose(kapa1);
     53    if(kapa2 != -1) KiiClose(kapa2);
    5254    return true;
    5355}
    … …  
    6062bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {
    6163    if (!pmVisualIsVisual() || !plotConvKernels) return true;
    62     if (!pmVisualInitWindow(&kapa, "ppSub:Images")) {
    63         return false;
    64     }
    65     pmVisualScaleImage(kapa, convKernels, "Convolution_Kernels", 0, true);
     64    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) {
     65        return false;
     66    }
     67    pmVisualScaleImage(kapa1, convKernels, "Convolution_Kernels", 0, true);
    6668    pmVisualAskUser(&plotConvKernels);
    6769    return true;
    … …  
    7476bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
    7577    if (!pmVisualIsVisual() || !plotStamps) return true;
    76     if (!pmVisualInitWindow (&kapa, "ppSub:Images")) {
     78    if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) {
    7779        return false;
    7880    }
    … …  
    134136        stampNum++;
    135137    }
    136     pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true);
     138    pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);
    137139
    138140    pmVisualAskUser(&plotStamps);
    … …  
    144146
    145147    if (!pmVisualIsVisual() || !plotLeastSquares) return true;
    146     if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
     148    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
    147149        return false;
    148150    }
    … …  
    157159        if (stamp == NULL) continue;
    158160
    159         psImage *im = stamp->matrix1;
    160         if (im == NULL) im = stamp->matrix2;
    161         if (im == NULL) im = stamp->matrixX;
     161        psImage *im = stamp->matrix;
    162162        if (im == NULL) continue;
    163163
    … …  
    187187        if (stamp == NULL) continue;
    188188
    189         psImage *im = stamp->matrix1;
    190         if (im == NULL) im = stamp->matrix2;
    191         if (im == NULL) im = stamp->matrixX;
     189        psImage *im = stamp->matrix;
    192190        if (im == NULL) continue;
    193191
    … …  
    197195
    198196    psImage *canvas32 = pmVisualImageToFloat(canvas);
    199     pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true);
     197    pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
    200198
    201199    pmVisualAskUser(&plotLeastSquares);
    … …  
    207205bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
    208206    if (!pmVisualIsVisual() || !plotImage) return true;
    209     if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
    210         return false;
    211     }
    212 
    213     pmVisualScaleImage(kapa, image, "Image", 0, true);
    214     pmVisualScaleImage(kapa, ref, "Reference", 1, true);
    215     pmVisualScaleImage(kapa, sub, "Subtraction", 2, true);
     207    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
     208        return false;
     209    }
     210
     211    pmVisualScaleImage(kapa1, image, "Image", 0, true);
     212    pmVisualScaleImage(kapa1, ref, "Reference", 1, true);
     213    pmVisualScaleImage(kapa1, sub, "Subtraction", 2, true);
     214    pmVisualAskUser(&plotImage);
     215    return true;
     216}
     217
     218bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) {
     219
     220    if (!pmVisualIsVisual()) return true;
     221    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
     222        return false;
     223    }
     224
     225    // get the kernel sizes
     226    int footprint = kernels->size;
     227
     228    // output image is a grid of NXsub by NYsub sub-images
     229    int NXsub = sqrt(kernels->num);
     230    int NYsub = kernels->num / NXsub;
     231    if (kernels->num % NXsub) NYsub++;
     232
     233    int NXpix = NXsub * (2*footprint + 1 + 3);
     234    int NYpix = NYsub * (2*footprint + 1 + 3);
     235
     236    psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
     237    psImageInit (output, 0.0);
     238
     239    for (int i = 0; i < kernels->num; i++) {
     240        pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i];
     241        psKernel *kernel = preCalc->kernel;
     242
     243        int xSub = i % NXsub;
     244        int ySub = i / NXsub;
     245
     246        int xPix = xSub * (2*footprint + 1 + 3) + footprint;
     247        int yPix = ySub * (2*footprint + 1 + 3) + footprint;
     248
     249        double sum = 0.0;
     250        for (int y = -footprint; y <= footprint; y++) {
     251            for (int x = -footprint; x <= footprint; x++) {
     252                output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
     253                sum += kernel->kernel[y][x];
     254            }
     255        }
     256        fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     257    }                                                   
     258       
     259    pmVisualScaleImage(kapa1, output, "Image", 0, true);
     260    pmVisualAskUser(&plotImage);
     261    return true;
     262}
     263
     264bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps) {
     265
     266    if (!pmVisualIsVisual()) return true;
     267    if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) {
     268        return false;
     269    }
     270
     271    // get the kernel sizes
     272    int footprint = stamps->footprint;
     273
     274    // choose the brightest stamp
     275    pmSubtractionStamp *maxStamp = NULL;
     276    float maxFlux = NAN;
     277    for (int i = 0; i < stamps->num; i++) {
     278        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     279        if (!isfinite(stamp->flux)) continue;
     280        if (!stamp->convolutions1 && !stamp->convolutions2) continue;
     281        if (!maxStamp) {
     282            maxFlux = stamp->flux;
     283            maxStamp = stamp;
     284            continue;
     285        }
     286        if (stamp->flux > maxFlux) {
     287            maxFlux = stamp->flux;
     288            maxStamp = stamp;
     289        }
     290    }
     291
     292    if (!isfinite(maxStamp->flux)) {
     293        fprintf (stderr, "no valid stamps?\n");
     294    }
     295
     296    int nKernels = 0;
     297
     298    if (maxStamp->convolutions1) {
     299        // output image is a grid of NXsub by NYsub sub-images
     300        nKernels = maxStamp->convolutions1->n;
     301        int NXsub = sqrt(nKernels);
     302        int NYsub = nKernels / NXsub;
     303        if (nKernels % NXsub) NYsub++;
     304
     305        int NXpix = NXsub * (2*footprint + 1 + 3);
     306        int NYpix = NYsub * (2*footprint + 1 + 3);
     307
     308        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
     309        psImageInit (output, 0.0);
     310
     311        for (int i = 0; i < nKernels; i++) {
     312            psKernel *kernel = maxStamp->convolutions1->data[i];
     313           
     314            int xSub = i % NXsub;
     315            int ySub = i / NXsub;
     316           
     317            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
     318            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
     319           
     320            double sum = 0.0;
     321            for (int y = -footprint; y <= footprint; y++) {
     322                for (int x = -footprint; x <= footprint; x++) {
     323                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
     324                    sum += kernel->kernel[y][x];
     325                }
     326            }
     327            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     328        }               
     329        pmVisualScaleImage(kapa2, output, "Image", 0, true);
     330    }                                   
     331       
     332    if (maxStamp->convolutions2) {
     333        // output image is a grid of NXsub by NYsub sub-images
     334        nKernels = maxStamp->convolutions2->n;
     335        int NXsub = sqrt(nKernels);
     336        int NYsub = nKernels / NXsub;
     337        if (nKernels % NXsub) NYsub++;
     338
     339        int NXpix = NXsub * (2*footprint + 1 + 3);
     340        int NYpix = NYsub * (2*footprint + 1 + 3);
     341
     342        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
     343        psImageInit (output, 0.0);
     344
     345        for (int i = 0; i < nKernels; i++) {
     346            psKernel *kernel = maxStamp->convolutions2->data[i];
     347           
     348            int xSub = i % NXsub;
     349            int ySub = i / NXsub;
     350           
     351            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
     352            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
     353           
     354            double sum = 0.0;
     355            for (int y = -footprint; y <= footprint; y++) {
     356                for (int x = -footprint; x <= footprint; x++) {
     357                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
     358                    sum += kernel->kernel[y][x];
     359                }
     360            }
     361            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     362        }               
     363        pmVisualScaleImage(kapa2, output, "Image", 1, true);
     364    }                                   
     365       
    216366    pmVisualAskUser(&plotImage);
    217367    return true;
    … …  
    240390
    241391        overlay[Noverlay].type = KII_OVERLAY_BOX;
     392        if ((stamp->x < 1.0) && (stamp->y < 1.0)) {
     393            // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
     394            continue;
     395        }
     396        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
     397            // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
     398            continue;
     399        }
    242400        overlay[Noverlay].x = stamp->x;
    243401        overlay[Noverlay].y = stamp->y;
    … …  
    255413}
    256414
     415static int footprint = 0;
     416static int NX = 0;
     417static int NY = 0;
     418static psImage *sourceImage      = NULL;
     419static psImage *targetImage      = NULL;
     420static psImage *residualImage    = NULL;
     421static psImage *fresidualImage   = NULL;
     422static psImage *differenceImage  = NULL;
     423static psImage *convolutionImage = NULL;
     424
     425bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
     426
     427    if (!pmVisualIsVisual()) return true;
     428
     429    // generate 4 storage images large enough to hold the N stamps:
     430
     431    footprint = stamps->footprint;
     432
     433    float NXf = sqrt(stamps->num);
     434    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
     435   
     436    float NYf = stamps->num / NX;
     437    NY = (int) NYf == NY ? NYf : NYf + 1.0;
     438
     439    int NXpix = (2*footprint + 1) * NX;
     440    NXpix += (NX > 1) ? 3 * NX : 0;
     441
     442    int NYpix = (2*footprint + 1) * NY;
     443    NYpix += (NY > 1) ? 3 * NY : 0;
     444
     445    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     446    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     447    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     448    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     449    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     450    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
     451   
     452    psImageInit (sourceImage,      0.0);
     453    psImageInit (targetImage,      0.0);
     454    psImageInit (residualImage,    0.0);
     455    psImageInit (fresidualImage,   0.0);
     456    psImageInit (differenceImage,  0.0);
     457    psImageInit (convolutionImage, 0.0);
     458
     459    return true;
     460}
     461
     462bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
     463
     464    if (!pmVisualIsVisual()) return true;
     465
     466    double sum;
     467
     468    int NXoff = index % NX;
     469    int NYoff = index / NX;
     470
     471    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
     472    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
     473
     474    // insert the (target) kernel into the (target) image:
     475    sum = 0.0;
     476    for (int y = -footprint; y <= footprint; y++) {
     477        for (int x = -footprint; x <= footprint; x++) {
     478            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
     479            sum += targetImage->data.F32[y + NYpix][x + NXpix];
     480        }
     481    }
     482    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     483
     484    // insert the (source) kernel into the (source) image:
     485    sum = 0.0;
     486    for (int y = -footprint; y <= footprint; y++) {
     487        for (int x = -footprint; x <= footprint; x++) {
     488            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
     489            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
     490        }
     491    }
     492    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     493
     494    // insert the (convolution) kernel into the (convolution) image:
     495    sum = 0.0;
     496    for (int y = -footprint; y <= footprint; y++) {
     497        for (int x = -footprint; x <= footprint; x++) {
     498            convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x];
     499            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
     500        }
     501    }
     502    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     503   
     504    // insert the (difference) kernel into the (difference) image:
     505    sum = 0.0;
     506    for (int y = -footprint; y <= footprint; y++) {
     507        for (int x = -footprint; x <= footprint; x++) {
     508            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
     509            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
     510        }
     511    }
     512    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     513
     514    // insert the (residual) kernel into the (residual) image:
     515    sum = 0.0;
     516    for (int y = -footprint; y <= footprint; y++) {
     517        for (int x = -footprint; x <= footprint; x++) {
     518            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
     519            sum += residualImage->data.F32[y + NYpix][x + NXpix];
     520        }
     521    }
     522    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     523
     524    // insert the (fresidual) kernel into the (fresidual) image:
     525    for (int y = -footprint; y <= footprint; y++) {
     526        for (int x = -footprint; x <= footprint; x++) {
     527            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
     528        }
     529    }
     530    return true;
     531}
     532
     533bool pmSubtractionVisualShowFit(double norm) {
     534
     535    // for testing, dump the residual image and exit
     536    if (0) {
     537        psMetadata *header = psMetadataAlloc();
     538        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
     539        psFits *fits = psFitsOpen("resid.fits", "w");
     540        psFitsWriteImage(fits, header, residualImage, 0, NULL);
     541        psFitsClose(fits);
     542        // exit (0);
     543    }
     544
     545    if (!pmVisualIsVisual()) return true;
     546    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
     547    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
     548
     549    KiiEraseOverlay (kapa1, "red");
     550    KiiEraseOverlay (kapa2, "red");
     551
     552    pmVisualScaleImage(kapa1, targetImage, "Target Stamps", 0, true);
     553    pmVisualScaleImage(kapa1, sourceImage, "Source Stamps", 1, true);
     554    pmVisualScaleImage(kapa1, convolutionImage, "Convolution Stamps", 2, true);
     555    KiiCenter (kapa1, 0.5*targetImage->numCols, 0.5*targetImage->numRows, 1);
     556
     557    pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true);
     558    pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);
     559
     560    if (1) {
     561        KiiImage image;
     562        KapaImageData data;
     563        Coords coords;
     564        strcpy (coords.ctype, "RA---TAN");
     565
     566        image.data2d = residualImage->data.F32;
     567        image.Nx = residualImage->numCols;
     568        image.Ny = residualImage->numRows;
     569        strcpy (data.name, "Residual Stamps");
     570        strcpy (data.file, "Residual Stamps");
     571
     572        data.zero  = -300.0;
     573        data.range = +600.0;
     574        data.logflux = 0;
     575
     576        KiiSetChannel (kapa2, 1);
     577        KiiNewPicture2D (kapa2, &image, &data, &coords);
     578    } else {
     579        pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true);
     580    }
     581
     582    KiiCenter (kapa2, 0.5*residualImage->numCols, 0.5*residualImage->numRows, 1);
     583
     584    pmVisualAskUser(NULL);
     585
     586    psFree(targetImage);
     587    psFree(sourceImage);
     588    psFree(convolutionImage);
     589    psFree(differenceImage);
     590    psFree(residualImage);
     591    psFree(fresidualImage);
     592
     593    targetImage = NULL;
     594    sourceImage = NULL;
     595    convolutionImage = NULL;
     596    differenceImage = NULL;
     597    residualImage = NULL;
     598    fresidualImage = NULL;
     599
     600    return true;
     601}
     602
     603bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels) {
     604
     605    Graphdata graphdata;
     606
     607    if (!pmVisualIsVisual()) return true;
     608    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;
     609
     610    KapaClearSections (kapa3);
     611    KapaInitGraph (&graphdata);
     612
     613    psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
     614    psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
     615
     616    graphdata.xmin = -1.0;
     617    graphdata.xmax = kernels->num + 1.0;
     618    graphdata.ymin = +32.0;
     619    graphdata.ymax = -32.0;
     620
     621    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0);
     622
     623    // construct the plot vectors
     624    for (int i = 0; i < kernels->num; i++) {
     625        x->data.F32[i] = i;
     626        y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
     627        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
     628        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
     629    }
     630    x->n = y->n = kernels->num;
     631
     632    float range;
     633    range = graphdata.xmax - graphdata.xmin;
     634    graphdata.xmax += 0.05*range;
     635    graphdata.xmin -= 0.05*range;
     636    range = graphdata.ymax - graphdata.ymin;
     637    graphdata.ymax += 0.05*range;
     638    graphdata.ymin -= 0.05*range;
     639
     640    KapaSetLimits (kapa3, &graphdata);
     641
     642    KapaSetFont (kapa3, "helvetica", 14);
     643    KapaBox (kapa3, &graphdata);
     644    KapaSendLabel (kapa3, "kernel number", KAPA_LABEL_XM);
     645    KapaSendLabel (kapa3, "coeff", KAPA_LABEL_YM);
     646
     647    graphdata.color = KapaColorByName ("black");
     648    graphdata.ptype = 2;
     649    graphdata.size = 0.5;
     650    graphdata.style = 2;
     651
     652    KapaPrepPlot   (kapa3, x->n, &graphdata);
     653    KapaPlotVector (kapa3, x->n, x->data.F32, "x");
     654    KapaPlotVector (kapa3, x->n, y->data.F32, "y");
     655
     656    psFree (x);
     657    psFree (y);
     658
     659    pmVisualAskUser(NULL);
     660    return true;
     661}
     662
    257663#else
    258664bool pmSubtractionVisualClose(void) {return true;}
    … …  
    261667bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps) {return true;}
    262668bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {return true;}
     669bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {return true;}
     670bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {return true;}
     671bool pmSubtractionVisualShowFit() {return true;}
     672bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
    263673#endif
Note: See TracChangeset for help on using the changeset viewer.