IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20335


Ignore:
Timestamp:
Oct 22, 2008, 1:35:13 PM (18 years ago)
Author:
beaumont
Message:

Added mosaic visualizations

Location:
branches/cnb_branch_20081011
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branch_20081011/ippconfig/recipes/psastro.config

    r19499 r20335  
    7575PSASTRO.MAX.INST.MAG.RAW       F32      0.0   # max instrumental magnitude for stars accepted for fitting
    7676
    77 PSASTRO.MATCH.LUMFUNC  BOOL     FALSE
     77PSASTRO.MATCH.LUMFUNC  BOOL     TRUE
    7878
    7979# option may be MAX, MIN, or VALUE. for VALUE, look up

  • branches/cnb_branch_20081011/psModules/src/astrom/Makefile.am

    r17036 r20335  
    1010        pmAstrometryModel.c \
    1111        pmAstrometryRefstars.c \
    12         pmAstrometryWCS.c
     12        pmAstrometryWCS.c \
     13        pmAstrometryVisual.c
    1314
    1415pkginclude_HEADERS = \
    … …  
    1920        pmAstrometryModel.h \
    2021        pmAstrometryRefstars.h \
    21         pmAstrometryWCS.h
     22        pmAstrometryWCS.h \
     23        pmAstrometryVisual.h
    2224
    2325CLEANFILES = *~

  • branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.c

    r20039 r20335  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2008-10-10 02:33:35 $
     10*  @version $Revision: 1.41.4.1 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2008-10-22 23:35:13 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2929#include <unistd.h>   // for unlink
    3030#include <pslib.h>
    31 #include <kapa.h> // cnb: do I need an IFDEF?
    3231
    3332#include "pmHDU.h"
    … …  
    3534#include "pmAstrometryObjects.h"
    3635#include "pmKapaPlots.h"
     36#include "pmAstrometryVisual.h"
    3737
    3838#define PM_ASTROMETRYOBJECTS_DEBUG 1
    … …  
    564564    psF32 **D2 = gridD2->data.F32;
    565565
    566 #ifdef PLOTS
    567     // vectors to hold dX and dY
    568     int nplot = raw->n * ref->n;
    569     float dXplot[nplot];
    570     float dYplot[nplot];
    571 #endif
    572566
    573567    // accumulate grids for focal plane (L,M) matches
    … …  
    579573            dY = ob1->FP->y - ob2->FP->y;
    580574
    581 #ifdef PLOTS
    582             dXplot[(i * ref->n) + j] = dX;
    583             dYplot[(i * ref->n) + j] = dY;
    584             // fprintf (f, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP->x, ob2->FP->x, ob1->FP->y, ob2->FP->y);
    585 #endif
    586575            // find bin coordinates for this delta-delta
    587576            if (!AstromGridBin (&iX, &iY, dX, dY)) {
    … …  
    618607        }
    619608
    620         # if 0
     609# if 0
    621610        char line[16];
    622611        psFits *fits = psFitsOpen ("grid.image.fits", "w");
    … …  
    625614        fprintf (stderr, "wrote grid image, press return to continue\n");
    626615        fgets (line, 15, stdin);
    627         # endif
     616# endif
    628617
    629618        // only check bins with at least 1/2 of max bin
    … …  
    670659        }
    671660
     661        pmAstromVisualPlotGridMatch (raw, ref,gridNP, stats->offset.x, stats->offset.y,
     662                                     maxOffpix, Scale, Offset);
     663
    672664        psFree (imStats);
    673665        // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better
    … …  
    689681    // fprintf (stderr, "sigma: nMatch: %d, nTest: %d, nTen: %d\n", stats->nMatch, stats->nTest, sort->data.U32[sort->n - 10]);
    690682
    691 #ifdef PLOTS
    692     /*** plot DXs and DYs ***/
    693     KapaSection section = {"s1", 0.00, 0.00, 1.0, 1.0};
    694     //KapaSection left = {"s2", 0.05, 0.05, 0.35, 0.8};
    695     //KapaSection right = {"s3", 0.55, 0.05, 0.35, 0.8};
    696     //KapaSection yprof = {"s4", 0.45, 0.05, 0.10, 0.9};
    697     //KapaSection xprof = {"s5", .05, .90, .35, .1};
    698 
    699     Graphdata graphdata;
    700     int kapa = pmKapaOpen(true);
    701     KapaClearPlots(kapa);
    702     KapaInitGraph(&graphdata);
    703     KapaSetSection(kapa, &section);
    704     graphdata.xmin = stats->offset.x - 1.5 * maxOffpix;
    705     graphdata.xmax = stats->offset.x + 1.5 * maxOffpix;
    706     graphdata.ymin = stats->offset.y - 1.5 * maxOffpix;
    707     graphdata.ymax = stats->offset.y + 1.5 * maxOffpix;
    708 
    709     KapaSetLimits(kapa, &graphdata);
    710     KapaSetFont(kapa, "helvetica", 14);
    711     KapaBox(kapa, &graphdata);
    712     KapaSendLabel (kapa, "X offset", KAPA_LABEL_XM);
    713     KapaSendLabel (kapa, "Y offset", KAPA_LABEL_YM);
    714     KapaSendLabel (kapa, "pmAstromGridAngle residuals. Big Box: Serach Region. Small Box: Correlation Peak.",
    715                    KAPA_LABEL_XP);
    716     graphdata.style = 2;
    717     graphdata.ptype = 0;
    718     graphdata.size = 0.4;
    719     graphdata.color = KapaColorByName ("black");
    720     KapaPrepPlot(kapa, nplot, &graphdata);
    721     KapaPlotVector (kapa, nplot, dXplot, "x");
    722     KapaPlotVector (kapa, nplot, dYplot, "y");
    723 
    724     //Overplot bounding box, peak of distribution
    725     float xbound[5] = { -maxOffpix, maxOffpix, maxOffpix, -maxOffpix, -maxOffpix};
    726     float ybound[5] = { -maxOffpix, -maxOffpix, maxOffpix, maxOffpix, -maxOffpix};
    727     float xbin[5] = {stats->offset.x - 0.5 * Scale, stats->offset.x + 0.5 * Scale, stats->offset.x + 0.5 * Scale,
    728                      stats->offset.x - 0.5 * Scale, stats->offset.x - 0.5 * Scale};
    729     float ybin[5] = {stats->offset.y - 0.5 * Scale, stats->offset.y - 0.5 * Scale, stats->offset.y + 0.5 * Scale,
    730                      stats->offset.y + 0.5 * Scale, stats->offset.y - 0.5 * Scale};
    731     graphdata.color = KapaColorByName("red");
    732     graphdata.style = 0;
    733     graphdata.size = 1.0;
    734     KapaPrepPlot(kapa, 5, &graphdata);
    735     KapaPlotVector (kapa, 5, xbound, "x");
    736     KapaPlotVector (kapa, 5, ybound, "y");
    737     KapaPrepPlot(kapa, 5, &graphdata);
    738     KapaPlotVector (kapa, 5, xbin, "x");
    739     KapaPlotVector (kapa, 5, ybin, "y");
    740 
    741     KapaPNG(kapa, "dXdY.png");
    742     char c;
    743     fscanf(stdin, "%c", &c);
    744     /*** end ***/
    745 #endif
    746683
    747684  psFree (sort);
    … …  
    767704    PS_ASSERT_PTR_NON_NULL(ref, NULL);
    768705    PS_ASSERT_PTR_NON_NULL(config, NULL);
    769 
    770706    bool status;
    771707    double xMin, xMax, yMin, yMax;

  • branches/cnb_branch_20081011/psModules/src/astrom/pmAstrometryObjects.h

    r15665 r20335  
    44 * @author EAM, IfA
    55 *
    6  * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-11-21 07:02:55 $
     6 * @version $Revision: 1.17.42.1 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2008-10-22 23:35:13 $
    88 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    99 */
    … …  
    238238 *
    239239 */
    240 bool pmAstromFitFPA(
    241     pmFPA *fpa,
    242     psArray *st1,
    243     psArray *st2,
    244     psArray *match,
    245     psMetadata *config
    246 );
     240bool pmAstromFitFPA(pmFPA *fpa,
     241                    psArray *st1,
     242                    psArray *st2,
     243                    psArray *match,
     244                    psMetadata *config);
    247245
    248246
    … …  
    269267    psArray *st2,
    270268    psArray *match,
    271     psMetadata *config
    272 );
     269    psMetadata *config);
    273270
    274271

  • branches/cnb_branch_20081011/psModules/src/psmodules.h

    r20067 r20335  
    8181#include <pmAstrometryRefstars.h>
    8282#include <pmAstrometryDistortion.h>
     83#include <pmAstrometryVisual.h>
    8384
    8485// the following headers are from psModule:imcombine

  • branches/cnb_branch_20081011/psastro/src/psastro.h

    r20263 r20335  
    5151psArray          *psastroRemoveClumps (psArray *input, int scale);
    5252
     53
    5354// utility functions:
    5455bool              psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars);
    … …  
    8485bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
    8586bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd);
     87bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
     88bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe);
    8689
    8790// demo plots

  • branches/cnb_branch_20081011/psastro/src/psastroMosaicOneChip.c

    r20067 r20335  
    55  if (!status) { \
    66   psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \
    7    return false; } 
     7   return false; }
    88
    99bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration) {
    … …  
    2929
    3030    // correct radius to FP units (physical pixel scale in microns per pixel)
    31     REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale"); 
     31    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale");
    3232
    3333    // allowed limits for valid solutions
    … …  
    4242    int order = psMetadataLookupS32 (&status, recipe, orderWord);
    4343    if (!status || (order == -1)) {
    44         order = defaultOrder;
     44        order = defaultOrder;
    4545    }
    4646
    … …  
    5050    if ((match->n <  8) && (order >= 1)) order = 0;
    5151    if ((match->n <  3) || (order < 0) || (order > 3)) {
    52         psLogMsg ("psastro", 3, "insufficient stars (%ld) or invalid order (%d)", match->n, order);
    53         return false;
    54     } 
     52        psLogMsg ("psastro", 3, "insufficient stars (%ld) or invalid order (%d)", match->n, order);
     53        return false;
     54    }
    5555
    5656    psLogMsg ("psastro", PS_LOG_DETAIL, "mosaic fit chip order %d", order);
    … …  
    6060    // coefficients frozen to the current values
    6161    if (order == 0) {
    62         // set FIT mask for all higher order terms of the existing solution
    63         // any existing SET masks will be retained.
    64         for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    65             for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    66                 if (i + j > 0) {
    67                     chip->toFPA->x->coeffMask[i][j] |= PS_POLY_MASK_FIT;
    68                     chip->toFPA->y->coeffMask[i][j] |= PS_POLY_MASK_FIT;
    69                 }
    70             }
    71         }
     62        // set FIT mask for all higher order terms of the existing solution
     63        // any existing SET masks will be retained.
     64        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     65            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     66                if (i + j > 0) {
     67                    chip->toFPA->x->coeffMask[i][j] |= PS_POLY_MASK_FIT;
     68                    chip->toFPA->y->coeffMask[i][j] |= PS_POLY_MASK_FIT;
     69                }
     70            }
     71        }
    7272    } else {
    73         psFree (chip->toFPA);
    74         chip->toFPA = psPlaneTransformAlloc (order, order);
    75         for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    76             for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    77                 if (i + j > order) {
    78                     chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
    79                     chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
    80                 }
    81             }
    82         }
     73        psFree (chip->toFPA);
     74        chip->toFPA = psPlaneTransformAlloc (order, order);
     75        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     76            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     77                if (i + j > order) {
     78                    chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
     79                    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
     80                }
     81            }
     82        }
    8383    }
    8484
    … …  
    8787    psStats *fitStats = NULL;
    8888//    if (order == 0) {
    89 //      fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    90 //      fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
    91 //      fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
     89//      fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     90//      fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
     91//      fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
    9292//    } else {
    93         fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    94         fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
    95         fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
     93        fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     94        fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
     95        fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
    9696//    }
    9797
    … …  
    9999    pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
    100100    if (!results) {
    101         psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n");
    102         return false;
     101        psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n");
     102        return false;
    103103    }
    104104
    … …  
    124124    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
    125125    if ((maxError > 0) && (astError > maxError)) {
    126         psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    127         validSolution = false;
     126        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
     127        validSolution = false;
    128128    }
    129129    if (astNstar < minNstar) {
    130         psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
    131         validSolution = false;
     130        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
     131        validSolution = false;
    132132    }
    133133
    … …  
    136136    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
    137137    if (validSolution) {
    138         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
    139         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     138        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
     139        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
    140140    } else {
    141         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
    142         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     141        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
     142        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
    143143    }
    144144    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
    … …  
    147147    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    148148
     149    //plot results
     150    psastroVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);
     151
    149152    psFree (fitStats);
    150153    psFree (results);

  • branches/cnb_branch_20081011/psastro/src/psastroMosaicSetMatch.c

    r17108 r20335  
    1010
    1111    // use small radius to match stars (assume starting astrometry is good)
    12     bool status = false; 
     12    bool status = false;
    1313    sprintf (radiusWord, "PSASTRO.MOSAIC.RADIUS.N%d", iteration);
    14     double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord); 
    15     if (!status) { 
    16         psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord);
    17         psFree (view);
    18         return false;
    19     } 
     14    double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord);
     15    if (!status) {
     16        psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord);
     17        psFree (view);
     18        return false;
     19    }
    2020
    2121    if (RADIUS <= 0.0) {
    22         if (iteration == 0) {
    23             psError(PS_ERR_IO, false, "Invalid match radius for first iteration: %s", radiusWord);
    24             psFree (view);
    25             return false;
    26         }
    27         psWarning ("skipping match for iteration %d\n", iteration);
    28         psFree (view);
    29         return true;
     22        if (iteration == 0) {
     23            psError(PS_ERR_IO, false, "Invalid match radius for first iteration: %s", radiusWord);
     24            psFree (view);
     25            return false;
     26        }
     27        psWarning ("skipping match for iteration %d\n", iteration);
     28        psFree (view);
     29        return true;
    3030    }
    3131
    … …  
    3434        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    3535        if (!chip->process || !chip->file_exists) { continue; }
    36         if (!chip->fromFPA) { continue; }
    37        
    38         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
     36        if (!chip->fromFPA) { continue; }
     37
     38        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
    3939            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    4040            if (!cell->process || !cell->file_exists) { continue; }
    4141
    42             // process each of the readouts
    43             // XXX there can only be one readout per chip, right?
    44             while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
    45                 if (! readout->data_exists) { continue; }
     42            // process each of the readouts
     43            // XXX there can only be one readout per chip, right?
     44            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
     45                if (! readout->data_exists) { continue; }
    4646
    47                 // select the raw objects for this readout
    48                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    49                 if (rawstars == NULL) { continue; }
     47                // select the raw objects for this readout
     48                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     49                if (rawstars == NULL) { continue; }
    5050
    51                 // select the raw objects for this readout
    52                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    53                 if (refstars == NULL) { continue; }
    54                 psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
     51                // select the raw objects for this readout
     52                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     53                if (refstars == NULL) { continue; }
     54                psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
    5555
    56                 psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS);
    57                 psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
     56                psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS);
     57                psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
    5858
    59                 // XXX drop the old one
    60                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY | PS_META_REPLACE, "astrometry matches", matches);
    61                 psFree (matches);
    62             }
    63         }
     59                psastroVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);
     60
     61                // XXX drop the old one
     62                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY | PS_META_REPLACE, "astrometry matches", matches);
     63                psFree (matches);
     64            }
     65        }
    6466    }
    6567    psFree (view);

  • branches/cnb_branch_20081011/psastro/src/psastroVisual.c

    r20264 r20335  
    1 /***Diagnostic plots for psastro***/
    2 
     1/***********************************/
     2/***Diagnostic plots for psastro****/
     3/***********************************/
    34# include "psastroInternal.h"
     5# include <string.h>
     6# include <stdlib.h>
     7# include <stdio.h>
    48
    59# if (HAVE_KAPA)
    … …  
    1115//variables to determine when things are plotted
    1216static bool isVisual             = false;
    13 static bool plotRawStars         = false;
    14 static bool plotRefStars         = false;
    15 static bool plotLumFunc          = false;
     17static bool plotRawStars         = true;
     18static bool plotRefStars         = true;
     19static bool plotLumFunc          = true;
    1620static bool plotRemoveClumps     = true;
    17 static bool plotOneChipFit       = false;
    18 static bool plotAstromGuessCheck = false;
     21static bool plotOneChipFit       = true;
     22static bool plotAstromGuessCheck = true;
     23static bool plotMosaicMatches    = true;
     24static bool plotMosaicOneChip    = true;
    1925
    2026// variables to store plotting window indices
    21 static int kapa = -1; //for plots
    22 static int kapa2 = -1; //also for plots
     27static int kapa = -1;
     28static int kapa2 = -1;
    2329
    2430// helper routine prototypes
    … …  
    2733bool psastroVisualTripleOverplot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing);
    2834bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing);
    29 
    30 
     35bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title);
     36
     37
     38/*****************************/
    3139/** Initialization Routines **/
    32 
    33 
    34 /**
    35  *  start or stop plotting
    36  */
     40/*****************************/
     41
     42/**  start or stop plotting */
    3743bool psastroSetVisual (bool mode) {
    3844    isVisual = mode;
    … …  
    4147
    4248
    43 /**
    44  * open, name, and resize a window if necessary
    45  */
     49/** open, name, and resize a window if necessary*/
    4650bool psastroVisualInitWindow( int *kapid, char *name ) {
    4751    if (*kapid == -1) {
    … …  
    5862
    5963
    60 /**
    61  * ask the user how to proceed
    62  */
     64/** ask the user how to proceed*/
    6365bool psastroVisualAskUser( bool *plotflag ) {
    6466    char key[10];
    … …  
    7779
    7880
    79 /**
    80  * destroy windows at the end of a run
    81  */
     81/** destroy windows at the end of a run*/
    8282bool psastroVisualClose() {
    8383    if(kapa != -1)
    … …  
    8888}
    8989
     90
     91/***********************/
    9092/** Plotting Routines **/
    91 
    92 /*  psastroVisualPlotRawStars */
    93 //  Show raw star astrometry as determined from first pass fit
    94 //  Invoked from within psastroAstromGuess
     93/***********************/
     94
     95
     96/**  psastroVisualPlotRawStars
     97 * Plot raw stars as determined from first pass astrometry fit
     98 * Called within psastroAstromGeuss
     99 */
    95100bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe)
    96101{
    97 
     102    // make sure we want to plot this
    98103    if (!plotRawStars || !isVisual) return true;
     104
     105    //set up plot region
    99106    if (!psastroVisualInitWindow (&kapa, "psastro:plots"))
    100107        return false;
    101 
    102108    Graphdata graphdata;
    103109    KapaSection section;
    104 
    105     bool status = false;
    106     float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
    107     float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");
    108110
    109111    KapaInitGraph (&graphdata);
    … …  
    118120    section.dy = 0.4;
    119121
     122    //initialize and populate plotting vectors
     123    bool status = false;
     124    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
     125    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");
     126
    120127    psVector *xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
    121128    psVector *yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
    … …  
    128135    KapaSetSection (kapa, &section);
    129136    psFree (section.name);
     137
     138    //Chip coordinates
    130139    int n = 0;
    131140    for (int i = 0; i < rawstars->n; i++) {
    … …  
    147156    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
    148157
     158    //Focal Plane Coordinates
    149159    section.x = 0.5;
    150160    section.y = 0.5;
    … …  
    173183    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
    174184
     185    // Tangent Plane Coordinates
    175186    section.x = 0.0;
    176187    section.y = 0.0;
    … …  
    199210    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
    200211
     212    //sky coordinates
    201213    section.x = 0.5;
    202214    section.y = 0.0;
    … …  
    250262
    251263
    252 /* psastroVisualPlotRefStars */
    253 // plot the location of references stars over the entire fpa
    254 // invoked during psastroChooseRefStars
     264/** psastroVisualPlotRefStars
     265 * plot the location of references stars over the entire fpa
     266 * invoked during psastroChooseRefStars
     267 */
    255268bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe)
    256269{
     270    //make sure we want to plot this
    257271    if (!isVisual || !plotRefStars) return true;
    258272    Graphdata graphdata;
     273
     274    //set up plotting variables
    259275    if (!psastroVisualInitWindow (&kapa, "psastro:plots"))
    260276        return false;
    261 
    262     bool status = false;
    263     float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
    264     float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
    265277
    266278    KapaInitGraph (&graphdata);
    … …  
    271283    graphdata.size = 0.5;
    272284    graphdata.style = 2;
     285
     286    //initialize and populate plot vectors
     287    bool status = false;
     288    float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
     289    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
    273290
    274291    psVector *xVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
    … …  
    309326
    310327
    311 /* psastroVisualPlotLuminosityFunction */
    312 // Plot the two luminosity functions created.
    313 // Called from psastroRefStarSubset.
    314 // The luminosity functions are used to select a subset of reference stars,
    315 // so we plot the cutoff that defines this subset
     328/** psastroVisualPlotLuminosityFunction
     329 * Plot the two luminosity functions created within psastroRefStarSubset
     330 * The luminosity functions are used to select a subset of reference stars,
     331 * so we plot the cutoff that defines this subset
     332 */
    316333bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag,
    317334                                          pmLumFunc *lumFunc, pmLumFunc *rawFunc) {
    318335
     336    // make sure we want to plot this
    319337    if ( !isVisual || !plotLumFunc ) return true;
     338
     339    //set up plotting variables
    320340    if ( !psastroVisualInitWindow (&kapa, "psastro:plots"))
    321341        return false;
    322342
    323     //Set up graph window
    324343    Graphdata graphdata;
    325344    KapaSection section1 = {"s1", .1, .1, .85, .35};
    … …  
    396415
    397416
    398 /**
     417/** psastroVisualPlotRemoveClumps
    399418 * Plot the stars in a region, and indicate which stars are part of 'clumps'
    400419 * These stars are flagged during astrometric fitting, since dense regions are
    401420 * harder to cross-match than sparse ones. Called during psastroRemoveClumps.
    402421 */
    403 
    404422bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {
     423
     424    //make sure we want to plot this
    405425    if (!isVisual || !plotRemoveClumps) return true;
     426
     427    //set up plot variables
    406428    if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) return false;
    407429
    408430    KapaSection section;
    409431    Graphdata graphdata;
    410 
    411     //initialize graph information
    412432    section.x = 0;
    413433    section.dx = .9;
    … …  
    419439    KapaSetSection(kapa, &section);
    420440    psFree(section.name);
     441
    421442    graphdata.ptype = 7;
    422443    graphdata.size = 0.5;
    423444    graphdata.style = 2;
    424445    graphdata.color = KapaColorByName ("black");
    425 
     446    KapaClearPlots(kapa);
    426447
    427448    //set up plot vectors
    … …  
    486507        }
    487508    }
     509
     510    //ask for user input and finish
    488511    psastroVisualAskUser (&plotRemoveClumps);
    489512    psFree (xVec);
    … …  
    494517
    495518
    496 /* psastroVisualPlotOneChipFit */
    497 // assess the goodness of fit for a signle chip by
    498 // plotting the fit residuals
    499 // invoked during psastroOneChipFit
     519/** psastroVisualPlotOneChipFit
     520 * assess the goodness of fit for a signle chip by
     521 * plotting the fit residuals
     522 * invoked during psastroOneChipFit
     523 */
    500524bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {
     525
     526    //make sure we want to plot this
    501527    if (!isVisual || !plotOneChipFit)
    502528        return true;
    503529
    504   if (!psastroVisualInitWindow (&kapa, "psastro:plots"))
     530    //plot the residuals
     531    if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)"))
    505532        return false;
    506533
     534    //ask for user input and finish
     535    psastroVisualAskUser(&plotOneChipFit);
     536    return true;
     537}
     538
     539
     540/**
     541 *  Plots the fpa chip corners projected on to the tangential plane before and after
     542 *  the astrometry solution has been applied. In psastroAstromGuessCheck, the old corners
     543 *  are then fit to the new corners to get a sense at how far off the initial WCS info was
     544 *  in offset, rotation, and scale. This procedure also plots the residuals of the fit from
     545 *  old to new coordinates
     546 */
     547bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo,
     548                                        psVector *cornerPn, psVector *cornerQn,
     549                                        psVector *cornerPd, psVector *cornerQd) {
     550
     551    //make sure we want to plot this
     552    if (!isVisual || !plotAstromGuessCheck) return true;
     553
     554    //set up graph window
     555    if ( !psastroVisualInitWindow (&kapa, "psastro:plots"))
     556        return false;
     557    Graphdata graphdata;
     558    KapaSection section;
     559    KapaInitGraph(&graphdata);
     560    KapaClearPlots (kapa);
     561
     562    graphdata.color = KapaColorByName ("black");
     563    graphdata.ptype = 2;
     564    graphdata.size = 0.5;
     565    graphdata.style = 2;
     566
     567    section.dx = 0.4;
     568    section.dy = 0.4;
     569
     570    //Old Corners
     571    section.x = 0.30;
     572    section.y = 0.50;
     573    section.name = NULL;
     574    psStringAppend (&section.name, "a0");
     575    KapaSetSection (kapa, &section);
     576    psFree(section.name);
     577
     578    psastroVisualScaleGraphdata (&graphdata, cornerPo, cornerPo);
     579    KapaSetLimits (kapa, &graphdata);
     580    KapaBox (kapa, &graphdata);
     581    KapaSendLabel (kapa, "P (Pixels)", KAPA_LABEL_XM);
     582    KapaSendLabel (kapa, "Q (Pixels)", KAPA_LABEL_YM);
     583    KapaSendLabel (kapa,
     584                   "Fiducial Points in the Tangent Plane. Black: Initial Astrometry. Red: Final Astrometry",
     585                   KAPA_LABEL_XP);
     586    KapaPrepPlot (kapa, cornerPo->n, &graphdata);
     587    KapaPlotVector (kapa, cornerPo->n, cornerPo->data.F32, "x");
     588    KapaPlotVector (kapa, cornerQo->n, cornerQo->data.F32, "y");
     589
     590    // New Corners
     591    graphdata.color = KapaColorByName("red");
     592    graphdata.ptype = 7;
     593    graphdata.size = 1.5;
     594    KapaPrepPlot (kapa, cornerPn->n, &graphdata);
     595    KapaPlotVector (kapa, cornerPn->n, cornerPn->data.F32, "x");
     596    KapaPlotVector (kapa, cornerQn->n, cornerQn->data.F32, "y");
     597
     598    // Residuals
     599    psVector *xResid = psVectorAlloc(cornerPn->n, PS_DATA_F32);
     600    psVector *yResid = psVectorAlloc(cornerQn->n, PS_DATA_F32);
     601    for(int i=0; i < cornerPn->n; i++) {
     602        xResid->data.F32[i] = (cornerPd->data.F32[i]);
     603        yResid->data.F32[i] = (cornerQd->data.F32[i]);
     604    }
     605
     606    graphdata.color = KapaColorByName("black");
     607    graphdata.size=0.5;
     608    section.x = 0.3;
     609    section.y = 0.0;
     610    section.name = NULL;
     611    psStringAppend (&section.name, "a1");
     612    KapaSetSection (kapa, &section);
     613    psFree(section.name);
     614
     615    psastroVisualScaleGraphdata (&graphdata, xResid, yResid);
     616    KapaSetLimits (kapa, &graphdata);
     617    KapaBox (kapa, &graphdata);
     618    KapaSendLabel (kapa, "dP", KAPA_LABEL_XM);
     619    KapaSendLabel (kapa, "dQ", KAPA_LABEL_YM);
     620    KapaSendLabel (kapa,
     621                   "Residual of the Fit from the Initial Astrometry to the Final Astrometry",
     622                   KAPA_LABEL_XP);
     623    KapaPrepPlot (kapa, cornerPd->n, &graphdata);
     624    KapaPlotVector (kapa, cornerPd->n, xResid->data.F32, "x");
     625    KapaPlotVector (kapa, cornerQd->n, yResid->data.F32, "y");
     626
     627    psFree(xResid);
     628    psFree(yResid);
     629    psastroVisualAskUser (&plotAstromGuessCheck);
     630    return true;
     631}
     632
     633
     634/**
     635 * psastroVisualPlotMosaicOneChip (rawstars, refstars, match);
     636 *   plot the residuals between raw stars and ref stars after
     637 *   fitting in psastroMosaicOneChip
     638 */
     639bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars,
     640                                     psArray *match, psMetadata *recipe)
     641{
     642
     643    //make sure we want to plot this
     644    if (!isVisual || !plotMosaicOneChip) return false;
     645
     646    //plot the residuals
     647    if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode"))
     648        return false;
     649
     650    //ask for user input and finish
     651    psastroVisualAskUser(&plotMosaicOneChip);
     652
     653    return true;
     654}
     655
     656/** psastroVisualPlotMosaicMatches
     657 * Plot the matches between raw and reference stars during psastroVisualMosaicSetMatch
     658 */
     659bool psastroVisualPlotMosaicMatches( psArray *rawstars, psArray *refstars,
     660                                    psArray *match, int iteration,
     661                                    psMetadata *recipe)
     662{
     663    //make sure we want to plot this
     664    if (!isVisual || !plotMosaicMatches) return true;
     665    char *title;
     666    switch(iteration) {
     667      case 0:
     668        title = "Matches found during psastroMosaicSetMatch iteration 0";
     669        break;
     670      case 1 :
     671        title = "Matches found during psastroMosaicSetMatch iteration 1";
     672        break;
     673      case 2:
     674        title = "Matches found during psastroMosaicSetMatch iteration 2";
     675        break;
     676      case 3:
     677        title = "Matches found during psastroMosaicSetMatch iteration 3";
     678        break;
     679      case 4:
     680         title = "Matches found during psastroMosaicSetMatch iteration 4";
     681         break;
     682      default:
     683        title = "Matches found during psastroMosaicSetMatch";
     684        break;
     685    }
     686
     687    if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, title))
     688        return false;
     689
     690    //ask for user input
     691    psastroVisualAskUser (&plotMosaicMatches);
     692    return true;
     693}
     694
     695
     696/*********************/
     697/** Helper Routines **/
     698/*********************/
     699
     700
     701/** psastroVisualTriplePlot
     702 * plot 2 vectors whose point sizes are scaled by a third vector
     703 * autoscales the plot
     704 */
     705bool psastroVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing)
     706{
     707    psastroVisualScaleGraphdata (graphdata, xVec, yVec);
     708    //printf("%f %f %f %f \n",graphdata->xmin, graphdata->xmax, graphdata->ymin, graphdata->ymax);
     709    // set the scale vector
     710    psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32);
     711    psastroVisualCreateScaleVec (zVec, zScale, increasing);
     712
     713    KapaSetFont (kapid, "helvetica", 14);
     714    KapaSetLimits(kapid, graphdata);
     715    KapaBox (kapid, graphdata);
     716
     717    // the point size will be scaled from the z vector
     718    graphdata->size = -1;
     719    KapaPrepPlot (kapid, xVec->n, graphdata);
     720    KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x");
     721    KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y");
     722    KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z");
     723    psFree (zScale);
     724    return true;
     725}
     726
     727
     728/** psastroVisualTripleOverplot
     729 * same as above, but does not rescale the plot or redraw the bounding box
     730 */
     731bool psastroVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) {
     732
     733    // set the scale vector
     734    psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32);
     735    psastroVisualCreateScaleVec (zVec, zScale, increasing);
     736
     737    KapaSetFont (kapid, "helvetica", 14);
     738
     739    // the point size will be scaled from the z vector
     740    graphdata->size = -1;
     741    KapaPrepPlot (kapid, xVec->n, graphdata);
     742    KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x");
     743    KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y");
     744    KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z");
     745    psFree (zScale);
     746    return true;
     747}
     748
     749
     750/** psastroVisualCreateScaleVec
     751 * create a vector of plot point sizes from a vector
     752 */
     753bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing) {
     754    // set limits based on data values
     755    float zmin = +FLT_MAX;
     756    float zmax = -FLT_MAX;
     757
     758    for (int i = 0; i < zVec->n; i++) {
     759        zmin = PS_MIN (zmin, zVec->data.F32[i]);
     760        zmax = PS_MAX (zmax, zVec->data.F32[i]);
     761    }
     762
     763    float range = zmax - zmin;
     764    if (range == 0.0) {
     765        psVectorInit (zScale, 1.0);
     766    } else {
     767        for (int i = 0; i < zVec->n; i++) {
     768            if (increasing) {
     769                zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zVec->data.F32[i] - zmin)/range));
     770            } else {
     771                zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zmax - zVec->data.F32[i])/range));
     772            }
     773        }
     774    }
     775    return true;
     776}
     777
     778
     779/** psastroVisualScaleGraphdata
     780 * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to
     781 * prevent outliers from making te plot region too big.
     782 */
     783bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec) {
     784
     785    graphdata->xmin = +FLT_MAX;
     786    graphdata->xmax = -FLT_MAX;
     787    graphdata->ymin = +FLT_MAX;
     788    graphdata->ymax = -FLT_MAX;
     789
     790    //determine standard deviation of xVec and yVec
     791    psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     792    psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     793    psVectorStats (statsX, xVec, NULL, NULL, 0);
     794    psVectorStats (statsY, yVec, NULL, NULL, 0);
     795
     796    float xhi  = statsX->sampleMedian + 3 *statsX->sampleStdev;
     797    float xlo = statsX->sampleMedian - 3 *statsX->sampleStdev;
     798    float yhi = statsY->sampleMedian + 3 *statsY->sampleStdev;
     799    float ylo = statsY->sampleMedian - 3 *statsY->sampleStdev;
     800
     801    // abort if there is no good data
     802    if (!isfinite(xhi) || !isfinite(xlo) || !isfinite(yhi) || !isfinite(ylo)) {
     803        graphdata->xmin = -1;
     804        graphdata->ymin  = -1;
     805        graphdata->xmax = 1;
     806        graphdata->ymax = 1;
     807        psFree(statsX);
     808        psFree(statsY);
     809        return false;
     810    }
     811
     812    for(int i = 0; i < xVec->n; i++) {
     813        if (!isfinite(xVec->data.F32[i])) continue;
     814        if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue;
     815        graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
     816        graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
     817    }
     818
     819    for (int i = 0; i < yVec->n; i++) {
     820        if (!isfinite(xVec->data.F32[i])) continue;
     821        if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue;
     822        graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
     823        graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
     824    }
     825
     826    // add a whitespace border
     827    float range = graphdata->xmax - graphdata->xmin;
     828    if (range == 0) range = 1;
     829    graphdata->xmin -= .05 * range;
     830    graphdata->xmax += .05 * range;
     831
     832    range = graphdata->ymax - graphdata->ymin;
     833    if (range == 0) range = 1;
     834    graphdata->ymin -= .05 * range;
     835    graphdata->ymax += .05 * range;
     836
     837    psFree (statsX);
     838    psFree (statsY);
     839    return true;
     840}
     841
     842
     843/** psastroVisualResdiPlot
     844 * Plots the residuals between matched raw and reference stars. Used to assess the quality of an
     845 * astrometry solution
     846 */
     847bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title) {
     848
     849    //set up the first window
     850    if (!psastroVisualInitWindow( &kapa, "psastro:plots")) return false;
     851
     852    //initialize graph information
     853    Graphdata graphdata;
     854    KapaSection section;
     855
     856    KapaInitGraph (&graphdata);
     857    KapaClearPlots (kapa);
     858
     859    graphdata.color = KapaColorByName ("black");
     860    graphdata.ptype = 7;
     861    graphdata.size = 0.5;
     862    graphdata.style = 2;
     863
     864    section.dx = 0.4;
     865    section.dy = 0.4;
     866
     867    //initialize and populate the plotting vectors
    507868    bool status = false;
    508869    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
    … …  
    511872    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
    512873
    513     //initialize graph information
    514     Graphdata graphdata;
    515     KapaSection section;
    516 
    517     KapaInitGraph (&graphdata);
    518     KapaClearPlots (kapa);
    519 
    520     graphdata.color = KapaColorByName ("black");
    521     graphdata.ptype = 7;
    522     graphdata.size = 0.5;
    523     graphdata.style = 2;
    524 
    525     section.dx = 0.4;
    526     section.dy = 0.4;
    527 
    528874    psVector *xVec = psVectorAlloc (match->n, PS_TYPE_F32);
    529875    psVector *yVec = psVectorAlloc (match->n, PS_TYPE_F32);
    … …  
    6611007    psStringAppend (&section.name, "a5");
    6621008    KapaSetSection (kapa, &section);
    663     KapaSendLabel (kapa, "Single Chip Fit Residuals (Chip Coordinates)", KAPA_LABEL_XP);
     1009    KapaSendLabel (kapa, title, KAPA_LABEL_XP);
    6641010    psFree (section.name);
    6651011
    6661012
    667     // *** X vs Y plot (different window)
     1013    // X vs Y plot (different window)
    6681014    if (!psastroVisualInitWindow( &kapa2, "psastro:plots"))
    6691015        return false;
    … …  
    7741120    psastroVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);
    7751121
    776     psastroVisualAskUser (&plotOneChipFit);
    777 
    7781122    psFree (xVec);
    7791123    psFree (yVec);
    … …  
    7821126}
    7831127
    784 
    785 /**
    786  *  Plots the fpa chip corners projected on to the tangential plane before and after
    787  *  the astrometry solution has been applied. In psastroAstromGuessCheck, the old corners
    788  *  are then fit to the new corners to get a sense at how far off the initial WCS info was
    789  *  in offset, rotation, and scale. This procedure also plots the residuals of the fit from
    790  *  old to new coordinates
    791  */
    792 bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo,
    793                                         psVector *cornerPn, psVector *cornerQn,
    794                                         psVector *cornerPd, psVector *cornerQd) {
    795 
    796     if (!isVisual || !plotAstromGuessCheck) return true;
    797     if ( !psastroVisualInitWindow (&kapa, "psastro:plots"))
    798         return false;
    799 
    800     //set up graph window
    801     Graphdata graphdata;
    802     KapaSection section;
    803     KapaInitGraph(&graphdata);
    804     KapaClearPlots (kapa);
    805 
    806     graphdata.color = KapaColorByName ("black");
    807     graphdata.ptype = 2;
    808     graphdata.size = 0.5;
    809     graphdata.style = 2;
    810 
    811     section.dx = 0.4;
    812     section.dy = 0.4;
    813 
    814     //Old Corners
    815     section.x = 0.30;
    816     section.y = 0.50;
    817     section.name = NULL;
    818     psStringAppend (&section.name, "a0");
    819     KapaSetSection (kapa, &section);
    820     psFree(section.name);
    821 
    822     psastroVisualScaleGraphdata (&graphdata, cornerPo, cornerPo);
    823     KapaSetLimits (kapa, &graphdata);
    824     KapaBox (kapa, &graphdata);
    825     KapaSendLabel (kapa, "P (Pixels)", KAPA_LABEL_XM);
    826     KapaSendLabel (kapa, "Q (Pixels)", KAPA_LABEL_YM);
    827     KapaSendLabel (kapa,
    828                    "Fiducial Points in the Tangent Plane. Black: Initial Astrometry. Red: Final Astrometry",
    829                    KAPA_LABEL_XP);
    830     KapaPrepPlot (kapa, cornerPo->n, &graphdata);
    831     KapaPlotVector (kapa, cornerPo->n, cornerPo->data.F32, "x");
    832     KapaPlotVector (kapa, cornerQo->n, cornerQo->data.F32, "y");
    833 
    834     // New Corners
    835     graphdata.color = KapaColorByName("red");
    836     graphdata.ptype = 7;
    837     graphdata.size = 1.5;
    838     KapaPrepPlot (kapa, cornerPn->n, &graphdata);
    839     KapaPlotVector (kapa, cornerPn->n, cornerPn->data.F32, "x");
    840     KapaPlotVector (kapa, cornerQn->n, cornerQn->data.F32, "y");
    841 
    842     // Residuals
    843     psVector *xResid = psVectorAlloc(cornerPn->n, PS_DATA_F32);
    844     psVector *yResid = psVectorAlloc(cornerQn->n, PS_DATA_F32);
    845     for(int i=0; i < cornerPn->n; i++) {
    846         xResid->data.F32[i] = (cornerPd->data.F32[i] / cornerPn->data.F32[i]);
    847         yResid->data.F32[i] = (cornerQd->data.F32[i] / cornerQn->data.F32[i]);
    848     }
    849 
    850     graphdata.color = KapaColorByName("black");
    851     graphdata.size=0.5;
    852     section.x = 0.3;
    853     section.y = 0.0;
    854     section.name = NULL;
    855     psStringAppend (&section.name, "a1");
    856     KapaSetSection (kapa, &section);
    857     psFree(section.name);
    858 
    859     psastroVisualScaleGraphdata (&graphdata, xResid, yResid);
    860     KapaSetLimits (kapa, &graphdata);
    861     KapaBox (kapa, &graphdata);
    862     KapaSendLabel (kapa, "dP / P", KAPA_LABEL_XM);
    863     KapaSendLabel (kapa, "dQ / Q", KAPA_LABEL_YM);
    864     KapaSendLabel (kapa,
    865                    "Residual of the Fit from the Initial Astrometry to the Final Astrometry",
    866                    KAPA_LABEL_XP);
    867     KapaPrepPlot (kapa, cornerPd->n, &graphdata);
    868     KapaPlotVector (kapa, cornerPd->n, xResid->data.F32, "x");
    869     KapaPlotVector (kapa, cornerQd->n, yResid->data.F32, "y");
    870 
    871     psFree(xResid);
    872     psFree(yResid);
    873     psastroVisualAskUser (&plotAstromGuessCheck);
    874     return true;
    875 }
    876 
    877 /** Helper Routines **/
    878 
    879 
    880 /* psastroVisualTriplePlot */
    881 // plot 2 vectors whose point sizes are scaled by a third vector
    882 // autoscales the plot
    883 bool psastroVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing)
    884 {
    885     psastroVisualScaleGraphdata (graphdata, xVec, yVec);
    886 
    887     // set the scale vector
    888     psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32);
    889     psastroVisualCreateScaleVec (zVec, zScale, increasing);
    890 
    891     KapaSetFont (kapid, "helvetica", 14);
    892     KapaSetLimits(kapid, graphdata);
    893     KapaBox (kapid, graphdata);
    894 
    895     // the point size will be scaled from the z vector
    896     graphdata->size = -1;
    897     KapaPrepPlot (kapid, xVec->n, graphdata);
    898     KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x");
    899     KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y");
    900     KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z");
    901     psFree (zScale);
    902     return true;
    903 }
    904 
    905 
    906 /* psastroVisualTripleOverplot */
    907 // same as above, but does not rescale the plot or redraw the bounding box
    908 bool psastroVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) {
    909 
    910     // set the scale vector
    911     psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32);
    912     psastroVisualCreateScaleVec (zVec, zScale, increasing);
    913 
    914     KapaSetFont (kapid, "helvetica", 14);
    915 
    916     // the point size will be scaled from the z vector
    917     graphdata->size = -1;
    918     KapaPrepPlot (kapid, xVec->n, graphdata);
    919     KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x");
    920     KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y");
    921     KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z");
    922     psFree (zScale);
    923     return true;
    924 }
    925 
    926 
    927 /* psastroVisualCreateScaleVec */
    928 //create a vector of plot point sizes from a vector
    929 bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing) {
    930     // set limits based on data values
    931     float zmin = +FLT_MAX;
    932     float zmax = -FLT_MAX;
    933 
    934     for (int i = 0; i < zVec->n; i++) {
    935         zmin = PS_MIN (zmin, zVec->data.F32[i]);
    936         zmax = PS_MAX (zmax, zVec->data.F32[i]);
    937     }
    938 
    939     float range = zmax - zmin;
    940     if (range == 0.0) {
    941         psVectorInit (zScale, 1.0);
    942     } else {
    943         for (int i = 0; i < zVec->n; i++) {
    944             if (increasing) {
    945                 zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zVec->data.F32[i] - zmin)/range));
    946             } else {
    947                 zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zmax - zVec->data.F32[i])/range));
    948             }
    949         }
    950     }
    951     return true;
    952 }
    953 
    954 
    955 /** psastroVisualScaleGraphdata
    956  * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to
    957  * prevent outliers from making te plot region too big.
    958  */
    959 bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec) {
    960     graphdata->xmin = +FLT_MAX;
    961     graphdata->xmax = -FLT_MAX;
    962     graphdata->ymin = +FLT_MAX;
    963     graphdata->ymax = -FLT_MAX;
    964     bool baddataX = true;
    965     bool baddataY = true;
    966 
    967     //determine standard deviation of xVec and yVec
    968     psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    969     psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    970     psVectorStats (statsX, xVec, NULL, NULL, 0);
    971     psVectorStats (statsY, yVec, NULL, NULL, 0);
    972 
    973     float xhi  = statsX->sampleMedian + 3 *statsX->sampleStdev;
    974     float xlo = statsX->sampleMedian - 3 *statsX->sampleStdev;
    975     float yhi = statsY->sampleMedian + 3 *statsY->sampleStdev;
    976     float ylo = statsY->sampleMedian - 3 *statsY->sampleStdev;
    977 
    978     for(int i = 0; i < xVec->n; i++) {
    979         if (!isfinite(xVec->data.F32[i])) continue;
    980         if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue;
    981         baddataX = false;
    982         graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
    983         graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
    984     }
    985 
    986     for (int i = 0; i < yVec->n; i++) {
    987         if (!isfinite(yVec->data.F32[i])) continue;
    988         if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue;
    989         baddataY = false;
    990         graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
    991         graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
    992     }
    993 
    994     // check to see if there was no good data
    995     if ( baddataX || baddataY) {
    996         graphdata->xmin = -1;
    997         graphdata->ymin  = -1;
    998         graphdata->xmax = 1;
    999         graphdata->ymax = 1;
    1000         return false;
    1001     }
    1002     float range = graphdata->xmax - graphdata->xmin;
    1003     if (range == 0) range = 1;
    1004     graphdata->xmin -= .05 * range;
    1005     graphdata->xmax += .05 * range;
    1006 
    1007     range = graphdata->ymax - graphdata->ymin;
    1008     if (range == 0) range = 1;
    1009     graphdata->ymin -= .05 * range;
    1010     graphdata->ymax += .05 * range;
    1011 
    1012     psFree (statsX);
    1013     psFree (statsY);
    1014     return true;
    1015 }
    1016 
    10171128# endif
    10181129
Note: See TracChangeset for help on using the changeset viewer.