IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

Added mosaic visualizations

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.