IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21208


Ignore:
Timestamp:
Jan 28, 2009, 1:36:37 PM (17 years ago)
Author:
beaumont
Message:

Added QA plots to ppSub. Migrated psastro QA plots to psModules/extras.

Location:
branches/cnb_branch_20090113
Files:
4 added
1 deleted
25 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branch_20090113/ppSub/src/ppSub.c

    r21120 r21208  
    6363    psTimerStop();
    6464
     65    pmSubtractionVisualClose(); //close plot windows, if -visual is set
    6566    psFree(config);
    6667    pmModelClassCleanup();

  • branches/cnb_branch_20090113/ppSub/src/ppSubArguments.c

    r20568 r21208  
    244244    psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin2", 0, "Binning factor for second level", 0);
    245245    psMetadataAddStr(arguments, PS_LIST_TAIL, "-dumpconfig", 0, "file to dump configuration to", NULL);
     246    psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL);
    246247
    247248    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 4) {
    … …  
    357358    }
    358359
     360    if (psMetadataLookupBool(NULL, arguments, "-visual")) {
     361        pmSubtractionSetVisual(true);
     362    }
     363
    359364    // Translate the kernel type
    360365    psString type = psMetadataLookupStr(NULL, arguments, "-type"); // Name of kernel type

  • branches/cnb_branch_20090113/ppSub/src/ppSubReadout.c

    r21120 r21208  
    451451    outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask);
    452452
     453    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
     454
     455    // examine the subtraction
     456    pmSubtractionVisualShowSubtraction(minuend->image,
     457                                       subtrahend->image,
     458                                       outRO->image);
     459
    453460    psFree(inConv);
    454461    psFree(refConv);
    455 
    456     outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
    457462
    458463    // Copy astrometry over

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

    r20801 r21208  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2008-11-20 01:26:07 $
     10*  @version $Revision: 1.42.8.1 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2009-01-28 23:36:37 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3232#include "pmHDU.h"
    3333#include "pmFPA.h"
     34#include "pmFPAExtent.h"
     35#include "pmFPAfile.h"
    3436#include "pmAstrometryObjects.h"
    3537#include "pmKapaPlots.h"

  • branches/cnb_branch_20090113/psModules/src/astrom/pmAstrometryVisual.c

    r20801 r21208  
    1 /** Diagnostic plots called from within pmAstrometry routines.
    2  * @author      Chris Beaumont
    3  * @date        October 2008
    4  */
    5 
    6 /* Include Files   */
    71#ifdef HAVE_CONFIG_H
    82#include <config.h>
    … …  
    1610#include <pslib.h>
    1711
    18 #include "pmKapaPlots.h"
    19 
    2012#include "pmHDU.h"
    2113#include "pmFPA.h"
     14#include "pmFPAfile.h"
    2215#include "pmAstrometryObjects.h"
     16#include "pmAstrometryVisual.h"
     17#include "pmFPAExtent.h"
    2318
    2419# if (HAVE_KAPA)
    2520# include <kapa.h>
    26 
    27 # define KAPAX 700
    28 # define KAPAY 700
    29 
     21# include "pmKapaPlots.h"
     22#include "pmVisual.h"
    3023
    3124//variables to determine when things are plotted
    … …  
    3326static bool plotGridMatch        = true;
    3427static bool plotTweak            = true;
     28static bool plotRawStars         = true;
     29static bool plotRefStars         = false;
     30static bool plotLumFunc          = true;
     31static bool plotRemoveClumps     = true;
     32static bool plotOneChipFit       = true;
     33static bool plotFixChips         = true;
     34static bool plotAstromGuessCheck = true;
     35static bool plotMosaicMatches    = true;
     36static bool plotCommonScale      = true;
     37static bool plotMosaicOneChip    = true;
    3538
    3639// variables to store plotting window indices
    3740static int kapa = -1;
    38 
    39 
    40 /* Utility Routine Prototypes */
    41 bool pmAstromVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip);
     41static int kapa2 = -1;
    4242
    4343
    4444/* Initialization Routines  */
    4545
    46 
    47 /** start or stop plotting */
    4846bool pmAstromSetVisual (bool mode) {
    4947    isVisual = mode;
    … …  
    5250
    5351
    54 /** open name, and resize a window if necessary */
    55 bool pmAstromVisualInitWindow (int *kapid, char *name) {
    56     if (*kapid == -1) {
    57         *kapid = KapaOpenNamedSocket("kapa", name);
    58         if (*kapid == -1) {
    59             fprintf (stderr, "failure to open kapa; visual mode disabled for pmAstrom\n");
    60             isVisual = false;
    61             return false;
    62         }
    63         KapaResize (*kapid, KAPAX, KAPAY);
    64     }
    65     return true;
    66 }
    67 
    68 
    69 /** ask the user how to proceed */
    70 bool pmAstromVisualAskUser(bool *plotflag)
    71 {
    72     char key[10];
    73     fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?");
    74     if (!fgets(key, 8, stdin)) {
    75         psWarning("Unable to read option");
    76     }
    77     if (key[0] == 's') {
    78         *plotflag = false;
    79     }
    80     if (key[0] == 'a') {
    81         pmAstromSetVisual(false);
    82     }
    83     return true;
    84 }
    85 
    86 
    87 /** destroy windows at the end of a run*/
    8852bool pmAstromVisualClose()
    8953{
    9054    if(kapa != -1)
    9155        KiiClose(kapa);
     56    if(kapa2 != -1)
     57        KiiClose(kapa2);
    9258    return true;
    9359}
    9460
    9561
    96 /* plotting routines */
    97 
    98 /**
    99  * Plot the offset between every pair of reference and raw source locations. The peak of this
    100  * distribution nominally gives the offset, scale difference, and rotation of the two catalogs.
    101  * Overplots the location of this peak as determined by pmAstromGridAngle, as well as some profiles
    102  * along horizontal and vertical cuts through this peak.
    103  */
    104 bool pmAstromVisualPlotGridMatch (const psArray *raw, ///< raw stars
    105                                   const psArray *ref, ///< reference stars
    106                                   psImage *gridNP,    ///< a 2D histogram of raw-ref star distances
    107                                   double offsetX,     ///< The X location (FP coordinates) of the peak of gridNP
    108                                   double offsetY,     ///< the Y location (FP coordinates) of the peak of gridNP
    109                                   double maxOffpix,   ///< The half-width of gridNP in FP coordinates
    110                                   double Scale,       ///< The pixel size of gridNP in histogram-bin-coordinates
    111                                   double Offset       ///< The (x,y) location (histogram-bin coordinates) of the FP point (0,0) in gridNP
     62/* Plotting Routines */
     63bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe)
     64{
     65    // make sure we want to plot this
     66    if (!plotRawStars || !isVisual) return true;
     67
     68    //set up plot region
     69    if (!pmVisualInitWindow (&kapa, "psastro:plots"))
     70        return false;
     71    Graphdata graphdata;
     72    KapaSection section;
     73
     74    KapaInitGraph (&graphdata);
     75    KapaClearPlots (kapa);
     76
     77    graphdata.color = KapaColorByName ("black");
     78    graphdata.ptype = 7;
     79    graphdata.size = 0.5;
     80    graphdata.style = 2;
     81
     82    section.dx = 0.4;
     83    section.dy = 0.4;
     84
     85    //initialize and populate plotting vectors
     86    bool status = false;
     87    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
     88    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");
     89
     90    psVector *xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     91    psVector *yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     92    psVector *zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     93
     94    section.x = 0.0;
     95    section.y = 0.5;
     96    section.name = NULL;
     97    psStringAppend (&section.name, "a0");
     98    KapaSetSection (kapa, &section);
     99    psFree (section.name);
     100
     101    //Chip coordinates
     102    int n = 0;
     103    for (int i = 0; i < rawstars->n; i++) {
     104        pmAstromObj *raw = rawstars->data[i];
     105        if (!isfinite(raw->Mag)) continue;
     106        if (raw->Mag < iMagMin) continue;
     107        if (raw->Mag > iMagMax) continue;
     108
     109        xVec->data.F32[n] = raw->chip->x;
     110        yVec->data.F32[n] = raw->chip->y;
     111        zVec->data.F32[n] = raw->Mag;
     112        n++;
     113    }
     114    xVec->n = yVec->n = zVec->n = n;
     115
     116    KapaSendLabel(kapa, "Chip", KAPA_LABEL_XP);
     117    KapaSendLabel(kapa, "X", KAPA_LABEL_XM);
     118    KapaSendLabel(kapa, "Y", KAPA_LABEL_YM);
     119    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     120
     121    //Focal Plane Coordinates
     122    section.x = 0.5;
     123    section.y = 0.5;
     124    section.name = NULL;
     125    psStringAppend (&section.name, "a1");
     126    KapaSetSection (kapa, &section);
     127    psFree (section.name);
     128
     129    n = 0;
     130    for (int i = 0; i < rawstars->n; i++) {
     131        pmAstromObj *raw = rawstars->data[i];
     132        if (!isfinite(raw->Mag)) continue;
     133        if (raw->Mag < iMagMin) continue;
     134        if (raw->Mag > iMagMax) continue;
     135
     136        xVec->data.F32[n] = raw->FP->x;
     137        yVec->data.F32[n] = raw->FP->y;
     138        zVec->data.F32[n] = raw->Mag;
     139        n++;
     140    }
     141    xVec->n = yVec->n = zVec->n = n;
     142
     143    KapaSendLabel (kapa, "Focal Plane", KAPA_LABEL_XP);
     144    KapaSendLabel (kapa, "L", KAPA_LABEL_XM);
     145    KapaSendLabel (kapa, "M", KAPA_LABEL_YM);
     146    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     147
     148    // Tangent Plane Coordinates
     149    section.x = 0.0;
     150    section.y = 0.0;
     151    section.name = NULL;
     152    psStringAppend (&section.name, "a2");
     153    KapaSetSection (kapa, &section);
     154    psFree (section.name);
     155
     156    n = 0;
     157    for (int i = 0; i < rawstars->n; i++) {
     158        pmAstromObj *raw = rawstars->data[i];
     159        if (!isfinite(raw->Mag)) continue;
     160        if (raw->Mag < iMagMin) continue;
     161        if (raw->Mag > iMagMax) continue;
     162
     163        xVec->data.F32[n] = raw->TP->x;
     164        yVec->data.F32[n] = raw->TP->y;
     165        zVec->data.F32[n] = raw->Mag;
     166        n++;
     167    }
     168    xVec->n = yVec->n = zVec->n = n;
     169
     170    KapaSendLabel (kapa, "Tangential Plane", KAPA_LABEL_XP);
     171    KapaSendLabel (kapa, "P", KAPA_LABEL_XM);
     172    KapaSendLabel (kapa, "Q", KAPA_LABEL_YM);
     173    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     174
     175    //sky coordinates
     176    section.x = 0.5;
     177    section.y = 0.0;
     178    section.name = NULL;
     179    psStringAppend (&section.name, "a3");
     180    KapaSetSection (kapa, &section);
     181    psFree (section.name);
     182
     183    n = 0;
     184    for (int i = 0; i < rawstars->n; i++) {
     185        pmAstromObj *raw = rawstars->data[i];
     186        if (!isfinite(raw->Mag)) continue;
     187        if (raw->Mag < iMagMin) continue;
     188        if (raw->Mag > iMagMax) continue;
     189
     190        xVec->data.F32[n] = DEG_RAD*raw->sky->r;
     191        yVec->data.F32[n] = DEG_RAD*raw->sky->d;
     192        zVec->data.F32[n] = raw->Mag;
     193        n++;
     194    }
     195    xVec->n = yVec->n = zVec->n = n;
     196
     197    KapaSendLabel (kapa, "Sky", KAPA_LABEL_XP);
     198    KapaSendLabel(kapa, "RA", KAPA_LABEL_XM);
     199    KapaSendLabel(kapa, "Dec", KAPA_LABEL_YM);
     200    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     201
     202    // flip x (East increase to left)
     203    SWAP (graphdata.xmin, graphdata.xmax);
     204    KapaSetLimits (kapa, &graphdata);
     205
     206    // plot label
     207    section.x = 0.0;
     208    section.y = 0.0;
     209    section.dx = .95;
     210    section.dy = .95;
     211    section.name = NULL;
     212    psStringAppend (&section.name, "a5");
     213    KapaSetSection (kapa, &section);
     214    KapaSendLabel (kapa, "Raw Star Selection - Initial Astrometry", KAPA_LABEL_XP);
     215    psFree (section.name);
     216
     217    // pause and wait for user input:
     218    pmVisualAskUser(&plotRawStars, &isVisual);
     219
     220    psFree (xVec);
     221    psFree (yVec);
     222    psFree (zVec);
     223    return true;
     224}
     225
     226
     227bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe)
     228{
     229    //make sure we want to plot this
     230    if (!isVisual || !plotRefStars) return true;
     231
     232    //set up plotting variables
     233    if (!pmVisualInitWindow (&kapa, "psastro:plots"))
     234        return false;
     235
     236    Graphdata graphdata;
     237    KapaInitGraph (&graphdata);
     238    KapaClearSections (kapa);
     239
     240    graphdata.color = KapaColorByName ("black");
     241    graphdata.ptype = 7;
     242    graphdata.size = 0.5;
     243    graphdata.style = 2;
     244
     245    //initialize and populate plot vectors
     246    bool status = false;
     247    float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
     248    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
     249
     250    psVector *xVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     251    psVector *yVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     252    psVector *zVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     253
     254    int n = 0;
     255    for (int i = 0; i < refstars->n; i++) {
     256        pmAstromObj *ref = refstars->data[i];
     257        if (!isfinite(ref->Mag)) continue;
     258        if (ref->Mag > rMagMax) continue;
     259        if (ref->Mag < rMagMin) continue;
     260
     261        xVec->data.F32[n] = DEG_RAD*ref->sky->r;
     262        yVec->data.F32[n] = DEG_RAD*ref->sky->d;
     263        zVec->data.F32[n] = ref->Mag;
     264        n++;
     265    }
     266    xVec->n = yVec->n = zVec->n = n;
     267
     268    KapaSendLabel (kapa, "RA", KAPA_LABEL_XM);
     269    KapaSendLabel (kapa, "Dec", KAPA_LABEL_YM);
     270    KapaSendLabel (kapa, "Reference Stars", KAPA_LABEL_XP);
     271    pmVisualTriplePlot(kapa, &graphdata, xVec, yVec, zVec, false);
     272
     273    // flip x (East increase to left)
     274    SWAP (graphdata.xmin, graphdata.xmax);
     275    KapaSetLimits (kapa, &graphdata);
     276
     277    // pause and wait for user input:
     278    pmVisualAskUser(&plotRefStars, &isVisual);
     279
     280    psFree (xVec);
     281    psFree (yVec);
     282    psFree (zVec);
     283    return true;
     284}
     285
     286
     287bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag,   // Log(n) for each magnitude bin
     288                                          psVector *Mag,     // magnitude bins
     289                                          pmLumFunc *lumFunc,// Fit to the reference star luminosity function
     290                                          pmLumFunc *rawFunc // Fit to the raw star luminoisty function
     291                                          )
     292{
     293
     294    // make sure we want to plot this
     295    if ( !isVisual || !plotLumFunc ) return true;
     296
     297    //set up plotting variables
     298    if ( !pmVisualInitWindow (&kapa, "psastro:plots"))
     299        return false;
     300
     301    Graphdata graphdata;
     302    KapaSection section1 = {"s1", .1, .1, .85, .35};
     303    KapaSection section2 =  {"s2", .1, .5, .85, .35};
     304    KapaSection section;
     305    if(rawFunc == NULL)
     306        section = section1;
     307    else
     308        section = section2;
     309    if (rawFunc == NULL)
     310        KapaClearPlots(kapa);
     311    KapaInitGraph(&graphdata);
     312
     313    //Determine Plot Limits
     314    pmVisualScaleGraphdata(&graphdata, Mag, lnMag, false);
     315
     316    //Make a line for the fit
     317    float x[2] = {graphdata.xmin, graphdata.xmax};
     318    float y[2] = {lumFunc->offset + x[0] * lumFunc->slope,
     319                 lumFunc->offset + x[1] * lumFunc->slope};
     320
     321    //Plot Data
     322    KapaSetSection(kapa, &section);
     323    KapaSetLimits(kapa, &graphdata);
     324
     325    KapaSetFont (kapa, "helvetica", 14);
     326    KapaBox (kapa, &graphdata);
     327    KapaSendLabel (kapa, "Magnitude", KAPA_LABEL_XM);
     328    KapaSendLabel (kapa, "Log(N)", KAPA_LABEL_YM);
     329    if (rawFunc == NULL)
     330        KapaSendLabel (kapa, "Raw Star Luminosity Function", KAPA_LABEL_XP);
     331    else
     332        KapaSendLabel (kapa,
     333                       "Reference Star Luminosity Function, Shifted Raw Fit, and Cutoff",
     334                       KAPA_LABEL_XP);
     335    graphdata.color=KapaColorByName("black");
     336    graphdata.style = 1;
     337    KapaPrepPlot (kapa, lnMag->n, &graphdata);
     338    KapaPlotVector(kapa, lnMag->n,   Mag->data.F32, "x");
     339    KapaPlotVector(kapa, lnMag->n, lnMag->data.F32, "y");
     340
     341    //Overplot fit
     342    graphdata.style=0;
     343    KapaPrepPlot(kapa,2,&graphdata);
     344    KapaPlotVector(kapa, 2, x, "x");
     345    KapaPlotVector(kapa, 2, y, "y");
     346
     347    //If rawFunc was supplied, overplot the raw star fit + cutoff
     348    if( rawFunc != NULL) {
     349        double mRef = 0.5*(lumFunc->mMin + lumFunc->mMax);
     350        double logRho = mRef * lumFunc->slope + lumFunc->offset;
     351        double mRaw = (logRho - rawFunc->offset) / rawFunc->slope;
     352        double deltaM = mRef - mRaw;
     353        double mRefMax = rawFunc->mMax + deltaM;
     354
     355        float xraw[2] = {rawFunc->mMin + deltaM, rawFunc->mMax + deltaM};
     356        float yraw[2] = {rawFunc->offset + (rawFunc->slope) * rawFunc->mMin,
     357                        rawFunc->offset + (rawFunc->slope) * rawFunc->mMax};
     358        float x[2] = {mRefMax, mRefMax};
     359        float y[2] = {graphdata.ymin, graphdata.ymax};
     360        graphdata.color= KapaColorByName("red");
     361        KapaPrepPlot(kapa, 2, &graphdata);
     362        KapaPlotVector(kapa, 2, x, "x");
     363        KapaPlotVector(kapa, 2, y, "y");
     364        KapaPrepPlot (kapa, 2, &graphdata);
     365        KapaPlotVector (kapa, 2, xraw, "x");
     366        KapaPlotVector (kapa, 2, yraw, "y");
     367
     368        // pause and wait for user input:
     369        pmVisualAskUser (&plotLumFunc, &isVisual);
     370    }
     371    return true;
     372} // end of pmAstromVisualPlotLuminosityFunction
     373
     374
     375bool pmAstromVisualPlotRemoveClumps (psArray *input, // Array containing the field stars
     376                                    psImage *count, // A 2D histogram of the field star distribution
     377                                    int scale,      // The pixel size of the histogram
     378                                    float limit     // The minimum numuber of stars in a bin flagged as a clump
     379                                    )
     380{
     381
     382    //make sure we want to plot this
     383    if (!isVisual || !plotRemoveClumps) return true;
     384
     385    //set up plot variables
     386    if ( !pmVisualInitWindow (&kapa, "psastro:plots")) return false;
     387
     388    KapaSection section;
     389    Graphdata graphdata;
     390    section.x = 0;
     391    section.dx = .9;
     392    section.y = 0;
     393    section.dy = .9;
     394    section.name = NULL;
     395    psStringAppend( &section.name, "a0");
     396    KapaInitGraph(&graphdata);
     397    KapaSetSection(kapa, &section);
     398    psFree(section.name);
     399
     400    graphdata.ptype = 7;
     401    graphdata.size = 0.5;
     402    graphdata.style = 2;
     403    graphdata.color = KapaColorByName ("black");
     404    KapaClearPlots(kapa);
     405
     406    //set up plot vectors
     407    float Xmin = +FLT_MAX;
     408    float Xmax = -FLT_MAX;
     409    float Ymin = +FLT_MAX;
     410    float Ymax = -FLT_MAX;
     411    psVector *xVec = psVectorAlloc (input->n, PS_TYPE_F32);
     412    psVector *yVec = psVectorAlloc (input->n, PS_TYPE_F32);
     413
     414    //determine boundaries for histogram bin calculation
     415    int n = 0;
     416    for (int i=0; i< input->n; i++) {
     417        pmAstromObj *obj = (pmAstromObj *)input->data[i];
     418        if (!isfinite(obj->FP->x)) continue;
     419        if (!isfinite(obj->FP->y)) continue;
     420        xVec->data.F32[n] = obj->FP->x;
     421        yVec->data.F32[n] = obj->FP->y;
     422        Xmin = PS_MIN (Xmin, xVec->data.F32[n]);
     423        Xmax = PS_MAX (Xmax, xVec->data.F32[n]);
     424        Ymin = PS_MIN (Ymin, yVec->data.F32[n]);
     425        Ymax = PS_MAX (Ymax, yVec->data.F32[n]);
     426        n++;
     427    }
     428    xVec->n = yVec->n = n;
     429
     430    //plot stars
     431    graphdata.xmax = Xmax;
     432    graphdata.xmin = Xmin;
     433    graphdata.ymax = Ymax;
     434    graphdata.ymin = Ymin;
     435    KapaSetLimits (kapa, &graphdata);
     436    KapaSetFont (kapa, "helvetica", 14);
     437
     438    KapaBox (kapa, &graphdata);
     439
     440    KapaSendLabel (kapa, "L (pixels)", KAPA_LABEL_XM);
     441    KapaSendLabel (kapa, "M (pixels)", KAPA_LABEL_YM);
     442    KapaSendLabel (kapa, "Regions Flagged as Clumps (Red Boxes)",
     443                   KAPA_LABEL_XP);
     444
     445    KapaPrepPlot (kapa, xVec->n, &graphdata);
     446    KapaPlotVector (kapa, xVec->n, xVec->data.F32, "x");
     447    KapaPlotVector (kapa, yVec->n, yVec->data.F32, "y");
     448
     449    graphdata.color = KapaColorByName ("red");
     450    graphdata.style = 1;
     451
     452    //overplot clumpy regions excluded from analysis
     453    for(int i = 0; i < count->numCols; i++) {
     454        for (int j = 0; j < count->numRows; j++) {
     455            if(count->data.U32[j][i] <= limit) continue; //not a clump
     456            float Xbot = (i - 5) * scale + Xmin;
     457            float Ybot = (j - 5) * scale + Ymin;
     458            if(Xbot < graphdata.xmin || Xbot > graphdata.xmax ||
     459               Ybot < graphdata.ymin || Ybot > graphdata.ymax) continue;
     460            float x[5] = {Xbot, Xbot + scale, Xbot + scale, Xbot, Xbot};
     461            float y[5] = {Ybot, Ybot, Ybot + scale, Ybot + scale, Ybot};
     462            KapaPrepPlot (kapa, 5, &graphdata);
     463            KapaPlotVector (kapa, 5, x, "x");
     464            KapaPlotVector (kapa, 5, y, "y");
     465        }
     466    }
     467
     468    //ask for user input and finish
     469    pmVisualAskUser (&plotRemoveClumps, &isVisual);
     470    psFree (xVec);
     471    psFree (yVec);
     472
     473    return true;
     474}
     475
     476
     477bool pmAstromVisualPlotOneChipFit (psArray *rawstars, // stars detected in the image
     478                                  psArray *refstars, // reference stars over the same region
     479                                  psArray *match,    // contains which rawstars match to which refstars
     480                                  psMetadata *recipe // data reduction recipe
    112481                                  )
     482{
     483
     484    //make sure we want to plot this
     485    if (!isVisual || !plotOneChipFit)
     486        return true;
     487
     488    //plot the residuals
     489    if (!pmVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)", &kapa, &kapa2))
     490        return false;
     491
     492    //ask for user input and finish
     493    pmVisualAskUser(&plotOneChipFit, &isVisual);
     494    return true;
     495}
     496
     497
     498bool pmAstromVisualPlotFixChips (pmFPAfile *input, // focal plane array file
     499                                psVector *xOld, // old X location of chip cornerss
     500                                psVector *yOld // old Y location of chip corners
     501                                )
     502{
     503    //make sure we want to plot this
     504    if(!isVisual || !plotFixChips) return true;
     505
     506    if(!pmVisualInitWindow(&kapa, "psastro:plots")) return false;
     507
     508    KapaSection section = {"s1", .05, .05, .9, .9};
     509    Graphdata graphdata;
     510    KapaInitGraph( &graphdata);
     511    KapaClearPlots (kapa);
     512    graphdata.ptype = 2;
     513    graphdata.style = 2;
     514
     515    psVector *xNew = psVectorAllocEmpty (xOld->n, PS_TYPE_F32);
     516    psVector *yNew = psVectorAllocEmpty (yOld->n, PS_TYPE_F32);
     517
     518    // copy of the code in psastroFixChips that generated xOld, yOld, but for xNew, yNew
     519    pmFPAview *view = pmFPAviewAlloc (0);
     520
     521    pmChip *obsChip = NULL;
     522    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     523        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     524
     525        psRegion *region = pmChipPixels(obsChip);
     526        psPlane ptCP, ptFP;
     527
     528        ptCP.x = region->x0; ptCP.y = region->y0;
     529        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     530        psVectorAppend (xNew, ptFP.x);
     531        psVectorAppend (yNew, ptFP.y);
     532
     533        ptCP.x = region->x0; ptCP.y = region->y1;
     534        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     535        psVectorAppend (xNew, ptFP.x);
     536        psVectorAppend (yNew, ptFP.y);
     537
     538        ptCP.x = region->x1; ptCP.y = region->y1;
     539        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     540        psVectorAppend (xNew, ptFP.x);
     541        psVectorAppend (yNew, ptFP.y);
     542
     543        ptCP.x = region->x1; ptCP.y = region->y0;
     544        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     545        psVectorAppend (xNew, ptFP.x);
     546        psVectorAppend (yNew, ptFP.y);
     547
     548        psFree (region);
     549    }
     550
     551    //set up graph
     552    pmVisualScaleGraphdata(&graphdata, xOld, yOld, true);
     553    pmVisualInitGraph(kapa, &section, &graphdata);
     554    KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM);
     555    KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM);
     556    KapaSendLabel (kapa, "Chip corners before (black) and after (red) FixChips", KAPA_LABEL_XP);
     557    KapaPrepPlot (kapa, xOld->n, &graphdata);
     558    KapaPlotVector (kapa, xOld->n, xOld->data.F32, "x");
     559    KapaPlotVector (kapa, xOld->n, yOld->data.F32, "y");
     560
     561    graphdata.ptype = 1;
     562    graphdata.color = KapaColorByName("red");
     563    KapaPrepPlot (kapa, xNew->n, &graphdata);
     564    KapaPlotVector (kapa, xNew->n, xNew->data.F32, "x");
     565    KapaPlotVector (kapa, xNew->n, yNew->data.F32, "y");
     566
     567    pmVisualAskUser(&plotFixChips, &isVisual);
     568    psFree(xNew);
     569    psFree(yNew);
     570    psFree(view);
     571    return true;
     572}
     573
     574
     575bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, // P coordinates of chip corners before fitting
     576                                        psVector *cornerQo, // Q coordinates of chip corners before fitting
     577                                        psVector *cornerPn, // P coordinates of chip corners after fitting
     578                                        psVector *cornerQn, // Q coordinates of chip corners after fitting
     579                                        psVector *cornerPd, // P coordinate residuals of fit from old to new coordinates
     580                                        psVector *cornerQd  // Q coordinate residuals of fit from old to new coordinates
     581                                        )
     582{
     583
     584    //make sure we want to plot this
     585    if (!isVisual || !plotAstromGuessCheck) return true;
     586
     587    //set up graph window
     588    if ( !pmVisualInitWindow (&kapa, "psastro:plots"))
     589        return false;
     590    Graphdata graphdata;
     591    KapaSection section;
     592    KapaInitGraph(&graphdata);
     593    KapaClearPlots (kapa);
     594
     595    graphdata.color = KapaColorByName ("black");
     596    graphdata.ptype = 2;
     597    graphdata.size = 0.5;
     598    graphdata.style = 2;
     599
     600    section.dx = 0.4;
     601    section.dy = 0.4;
     602
     603    //Old Corners
     604    section.x = 0.30;
     605    section.y = 0.50;
     606    section.name = NULL;
     607    psStringAppend (&section.name, "a0");
     608    KapaSetSection (kapa, &section);
     609    psFree(section.name);
     610
     611    pmVisualScaleGraphdata (&graphdata, cornerPo, cornerPo, true);
     612    KapaSetLimits (kapa, &graphdata);
     613    KapaBox (kapa, &graphdata);
     614    KapaSendLabel (kapa, "P (Pixels)", KAPA_LABEL_XM);
     615    KapaSendLabel (kapa, "Q (Pixels)", KAPA_LABEL_YM);
     616    KapaSendLabel (kapa,
     617                   "Fiducial Points in the Tangent Plane. Black: Initial Astrometry. Red: Final Astrometry",
     618                   KAPA_LABEL_XP);
     619    KapaPrepPlot (kapa, cornerPo->n, &graphdata);
     620    KapaPlotVector (kapa, cornerPo->n, cornerPo->data.F32, "x");
     621    KapaPlotVector (kapa, cornerQo->n, cornerQo->data.F32, "y");
     622
     623    // New Corners
     624    graphdata.color = KapaColorByName("red");
     625    graphdata.ptype = 7;
     626    graphdata.size = 1.5;
     627    KapaPrepPlot (kapa, cornerPn->n, &graphdata);
     628    KapaPlotVector (kapa, cornerPn->n, cornerPn->data.F32, "x");
     629    KapaPlotVector (kapa, cornerQn->n, cornerQn->data.F32, "y");
     630
     631    // Residuals
     632    psVector *xResid = psVectorAlloc(cornerPn->n, PS_DATA_F32);
     633    psVector *yResid = psVectorAlloc(cornerQn->n, PS_DATA_F32);
     634    for(int i=0; i < cornerPn->n; i++) {
     635        xResid->data.F32[i] = (cornerPd->data.F32[i]);
     636        yResid->data.F32[i] = (cornerQd->data.F32[i]);
     637    }
     638
     639    graphdata.color = KapaColorByName("black");
     640    graphdata.size=0.5;
     641    section.x = 0.3;
     642    section.y = 0.0;
     643    section.name = NULL;
     644    psStringAppend (&section.name, "a1");
     645    KapaSetSection (kapa, &section);
     646    psFree(section.name);
     647
     648    pmVisualScaleGraphdata (&graphdata, xResid, yResid, true);
     649    KapaSetLimits (kapa, &graphdata);
     650    KapaBox (kapa, &graphdata);
     651    KapaSendLabel (kapa, "dP", KAPA_LABEL_XM);
     652    KapaSendLabel (kapa, "dQ", KAPA_LABEL_YM);
     653    KapaSendLabel (kapa,
     654                   "Residual of the Fit from the Initial Astrometry to the Final Astrometry",
     655                   KAPA_LABEL_XP);
     656    KapaPrepPlot (kapa, cornerPd->n, &graphdata);
     657    KapaPlotVector (kapa, cornerPd->n, xResid->data.F32, "x");
     658    KapaPlotVector (kapa, cornerQd->n, yResid->data.F32, "y");
     659
     660    psFree(xResid);
     661    psFree(yResid);
     662    pmVisualAskUser (&plotAstromGuessCheck, &isVisual);
     663    return true;
     664}
     665
     666
     667bool pmAstromVisualPlotCommonScale (pmFPA *fpa,         // the fpa
     668                                   psVector *oldScale  // the old pixel scale of each chip in the fpa
     669                                   )
     670{
     671    //make sure we want to plot this
     672    if (!isVisual || !plotCommonScale) return false;
     673
     674    if (!pmVisualInitWindow(&kapa, "psastro:plots")) return false;
     675
     676    KapaSection section = {"s1", .05, .05, .9, .9};
     677    Graphdata graphdata;
     678    psPlane ptCH, ptFP;
     679    ptCH.x = 0;
     680    ptCH.y = 0;
     681    psVector *xVec = psVectorAlloc (oldScale->n, PS_TYPE_F32);
     682    psVector *yVec = psVectorAlloc (oldScale->n, PS_TYPE_F32);
     683
     684    int nobj = 0;
     685
     686    //project each chip corner to the Focal Plane
     687    for(int i = 0; i < fpa->chips->n; i++) {
     688        pmChip *chip = fpa->chips->data[i];
     689        if (!chip->process || !chip->file_exists) { continue; }
     690        if (!chip->toFPA) { continue; }
     691
     692        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     693        xVec->data.F32[nobj] = ptFP.x;
     694        yVec->data.F32[nobj] = ptFP.y;
     695        nobj++;
     696    }
     697
     698    //set up plot window
     699    KapaInitGraph (&graphdata);
     700    KapaClearPlots (kapa);
     701    KapaSetSection (kapa, &section);
     702    KapaSetFont (kapa, "helvetica", 14);
     703    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, oldScale, false);
     704    KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM);
     705    KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM);
     706    KapaSendLabel (kapa, "Old Pixel Scale of FPA Chips (Not to Scale)", KAPA_LABEL_XP);
     707
     708    pmVisualAskUser (&plotCommonScale, &isVisual);
     709    return true;
     710}
     711
     712
     713bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars,
     714                                     psArray *match, psMetadata *recipe)
     715{
     716
     717    //make sure we want to plot this
     718    if (!isVisual || !plotMosaicOneChip) return false;
     719
     720    //plot the residuals
     721    if (!pmVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode", &kapa, &kapa2))
     722        return false;
     723
     724    //ask for user input and finish
     725    pmVisualAskUser(&plotMosaicOneChip, &isVisual);
     726
     727    return true;
     728}
     729
     730
     731bool pmAstromVisualPlotMosaicMatches( psArray *rawstars, psArray *refstars,
     732                                    psArray *match, int iteration,
     733                                    psMetadata *recipe)
     734{
     735    //make sure we want to plot this
     736    if (!isVisual || !plotMosaicMatches) return true;
     737
     738    char title[60];
     739    sprintf(title, "Matches found during psastroMosaicSetMatch iteration %d", iteration);
     740
     741    if (!pmVisualResidPlot(rawstars, refstars, match, recipe, title, &kapa, &kapa2))
     742        return false;
     743
     744    //ask for user input
     745    pmVisualAskUser (&plotMosaicMatches, &isVisual);
     746    return true;
     747}
     748
     749
     750bool pmAstromVisualPlotGridMatch (const psArray *raw,
     751                                  const psArray *ref,
     752                                  psImage *gridNP,
     753                                  double offsetX,
     754                                  double offsetY,
     755                                  double maxOffpix,
     756                                  double Scale,
     757                                  double Offset)
    113758{
    114759    //make sure we want to plot this
    115760    if (!isVisual || !plotGridMatch) return true;
    116     if (!pmAstromVisualInitWindow(&kapa, "pmAstrom:plots"))
     761    if (!pmVisualInitWindow(&kapa, "pmAstrom:plots"))
    117762        return false;
    118763
    … …  
    254899    KapaPlotVector (kapa, 2, yslice, "y");
    255900
    256     pmAstromVisualAskUser(&plotGridMatch);
     901    pmVisualAskUser(&plotGridMatch, &isVisual);
    257902    return true;
    258903} // end of pmAstromVisualPlotGridMatch
    259904
    260905
    261 /** Plot the refinements made within pmAstromGridTweak.
    262  * After pmAstromGridMatch finds the best rotaion/scale/offset between raw and reference stars
    263  * within a coarse grid of rotations/scales, pmAstromGridTweak computes a higher precision
    264  * estimate of the offset. It computes the 2 point correlation function between raw and ref
    265  * stars along horizontal and vertical cuts through the first-guess offset. It finds the peak
    266  * of these two profiles and adjusts the offset accordingly. This procedure plots the profiles.
    267  */
    268 bool pmAstromVisualPlotTweak (psVector *xHist, ///< Smoothed Horizontal cut through the histogram
    269                               psVector *yHist, ///< Smoothed Vertical cut throug the histogram
    270                               int xBin,        ///< X Bin index of the histogram peak
    271                               int yBin         ///< Y bin index of the histogram peak
     906bool pmAstromVisualPlotTweak (psVector *xHist, // Smoothed Horizontal cut through the histogram
     907                              psVector *yHist, // Smoothed Vertical cut throug the histogram
     908                              int xBin,        // X Bin index of the histogram peak
     909                              int yBin         // Y bin index of the histogram peak
    272910    )
    273911{
    274912    //make sure we want to plot this
    275913    if (!isVisual || !plotTweak) return true;
    276     if (!pmAstromVisualInitWindow(&kapa, "pmAstrom:plots")) {
     914    if (!pmVisualInitWindow(&kapa, "pmAstrom:plots")) {
    277915        return false;
    278916    }
    … …  
    298936
    299937    // plot the X histogram
    300     pmAstromVisualScaleGraphdata(&graphdata, xIndices, xHist, false);
     938    pmVisualScaleGraphdata(&graphdata, xIndices, xHist, false);
    301939    KapaSetSection(kapa, &section1);
    302940    KapaSetLimits (kapa, &graphdata);
    … …  
    325963
    326964    //plot the Y histogram
    327     pmAstromVisualScaleGraphdata(&graphdata, yIndices, yHist, false);
     965    pmVisualScaleGraphdata(&graphdata, yIndices, yHist, false);
    328966    KapaSetSection(kapa, &section2);
    329967    KapaSetLimits (kapa, &graphdata);
    … …  
    357995                   KAPA_LABEL_XP);
    358996
    359     pmAstromVisualAskUser(&plotTweak);
     997    pmVisualAskUser(&plotTweak, &isVisual);
    360998
    361999    psFree(xIndices);
    … …  
    3641002} //end of pmAstromPlotTweak
    3651003
    366 
    367 /** pmAstromScaleGraphdata
    368  * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to
    369  * prevent outliers from making te plot region too big.
    370  */
    371 bool pmAstromVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip) {
    372 
    373     graphdata->xmin = +FLT_MAX;
    374     graphdata->xmax = -FLT_MAX;
    375     graphdata->ymin = +FLT_MAX;
    376     graphdata->ymax = -FLT_MAX;
    377 
    378     //determine standard deviation of xVec and yVec
    379     psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    380     psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    381     psVectorStats (statsX, xVec, NULL, NULL, 0);
    382     psVectorStats (statsY, yVec, NULL, NULL, 0);
    383 
    384     float xhi  = statsX->sampleMedian + 3 *statsX->sampleStdev;
    385     float xlo = statsX->sampleMedian - 3 *statsX->sampleStdev;
    386     float yhi = statsY->sampleMedian + 3 *statsY->sampleStdev;
    387     float ylo = statsY->sampleMedian - 3 *statsY->sampleStdev;
    388 
    389     // don't sigma clip
    390     if (!clip) {
    391         xhi = +FLT_MAX;
    392         xlo = -FLT_MAX;
    393         yhi = +FLT_MAX;
    394         ylo = -FLT_MAX;
    395     }
    396 
    397     // abort if there is no good data
    398     if (!isfinite(xhi) || !isfinite(xlo) || !isfinite(yhi) || !isfinite(ylo)) {
    399         graphdata->xmin = -1;
    400         graphdata->ymin  = -1;
    401         graphdata->xmax = 1;
    402         graphdata->ymax = 1;
    403         psFree(statsX);
    404         psFree(statsY);
    405         return false;
    406     }
    407 
    408     for(int i = 0; i < xVec->n; i++) {
    409         if (!isfinite(xVec->data.F32[i])) continue;
    410         if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue;
    411         graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
    412         graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
    413     }
    414 
    415     for (int i = 0; i < yVec->n; i++) {
    416         if (!isfinite(xVec->data.F32[i])) continue;
    417         if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue;
    418         graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
    419         graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
    420     }
    421 
    422     // add a whitespace border
    423     float range = graphdata->xmax - graphdata->xmin;
    424     if (range == 0) range = 1;
    425     graphdata->xmin -= .05 * range;
    426     graphdata->xmax += .05 * range;
    427 
    428     range = graphdata->ymax - graphdata->ymin;
    429     if (range == 0) range = 1;
    430     graphdata->ymin -= .05 * range;
    431     graphdata->ymax += .05 * range;
    432 
    433     psFree (statsX);
    434     psFree (statsY);
    435     return true;
    436 }
    437 
    438 
    4391004# else
    4401005
    4411006bool pmAstromSetVisual(bool mode) { return true; }
    442 bool pmAstromVisualInitWindow (int *kapid, char *name) { return true; }
    443 bool pmAstromVisualAskUser (bool *plotflag) { return true; }
    4441007bool pmAstromVisualClose() { return true; }
    4451008bool pmAstromVisualPlotGridMatch (const psArray *raw, const psArray *ref, psImage *gridNP, double offsetX, double offsetY, double maxOffpix, double Scale, double Offset) { return true; }
    4461009bool pmAstromVisualPlotTweak (psVector *xHist, psVector *yHist, int xBin, int yBin) {return true;}
     1010bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc) {return true;}
     1011bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) {return true;}
     1012bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe) {return true;}
     1013bool pmAstromVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {return true;}
     1014bool pmAstromVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld) {return true;}
     1015bool pmAstromVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;}
     1016bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd) {return true;}
     1017bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;}
     1018bool pmAstromVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale) {return true;}
     1019bool pmAstromVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe) {return true;}
    4471020
    4481021# endif

  • branches/cnb_branch_20090113/psModules/src/astrom/pmAstrometryVisual.h

    r20801 r21208  
     1/*
     2 * @file pmAstrometryVisual.h
     3 * @author Chris Beaumont, IfA
     4 * @brief A set of functions to display visual diagnostics from psastro
     5 * Copyright 2009 Institute for Astronomy, University of Hawaii
     6 */
     7
     8#ifndef PM_ASTROM_VISUAL_H
     9#define PM_ASTROM_VISUAL_H
     10
     11
     12/** A fit to the logN / logS curve for a set of stars
     13 * logN = offset + slope * logS
     14 */
     15typedef struct {
     16    double mMin;                        ///< minimum magnitude bin with data
     17    double mMax;                        ///< maximum magnitude bin with data
     18    double offset;                      ///< fitted line offset
     19    double slope;                       ///< fitted line slope
     20    double mPeak;                       ///< mag of peak bin
     21    int nPeak;                          ///< # of stars in peak bin
     22    int sPeak;                          ///< sum of stars to peak bin
     23} pmLumFunc;
     24
     25
     26/** Enable or disable visual plotting for psastro routines
     27 * @param mode true/false to enable/disable plotting
     28 * @return true for success */
    129bool pmAstromSetVisual(bool mode);
    2 bool pmAstromVisualInitWindow (int *kapid, char *name);
    3 bool pmAstromVisualAskUser (bool *plotflag);
     30
     31
     32/** Close plotting windows at the end of a run
     33 * @return true for success */
    434bool pmAstromVisualClose();
    5 bool pmAstromVisualPlotGridMatch (const psArray *raw, const psArray *ref, psImage *gridNP, double offsetX, double offsetY, double maxOffpix, double Scale, double Offset);
    6 bool pmAstromVisualPlotTweak (psVector *xHist, psVector *yHist, int xBin, int yBin);
     35
     36
     37/**
     38 * Plot the offset between every pair of reference and raw source locations. The peak of this
     39 * distribution nominally gives the offset, scale difference, and rotation of the two catalogs.
     40 * Overplots the location of this peak as determined by pmAstromGridAngle, as well as some profiles
     41 * along horizontal and vertical cuts through this peak.
     42 */
     43bool pmAstromVisualPlotGridMatch (const psArray *raw, ///< raw stars
     44                                  const psArray *ref, ///< reference stars
     45                                  psImage *gridNP,    ///< a 2D histogram of raw-ref star distances
     46                                  double offsetX,     ///< The X location (FP coordinates) of the peak of gridNP
     47                                  double offsetY,     ///< the Y location (FP coordinates) of the peak of gridNP
     48                                  double maxOffpix,   ///< The half-width of gridNP in FP coordinates
     49                                  double Scale,       ///< The pixel size of gridNP in histogram-bin-coordinates
     50                                  double Offset       ///< The (x,y) location (histogram-bin coordinates) of the FP point (0,0) in gridNP
     51                                  );
     52
     53
     54/**
     55 * Plot the refinements made within pmAstromGridTweak.
     56 * After pmAstromGridMatch finds the best rotaion/scale/offset between raw and reference stars
     57 * within a coarse grid of rotations/scales, pmAstromGridTweak computes a higher precision
     58 * estimate of the offset. It computes the 2 point correlation function between raw and ref
     59 * stars along horizontal and vertical cuts through the first-guess offset. It finds the peak
     60 * of these two profiles and adjusts the offset accordingly. This procedure plots the profiles.
     61 */
     62bool pmAstromVisualPlotTweak (psVector *xHist, ///< Smoothed Horizontal cut through the histogram
     63                              psVector *yHist, ///< Smoothed Vertical cut throug the histogram
     64                              int xBin,        ///< X Bin index of the histogram peak
     65                              int yBin         ///< Y bin index of the histogram peak
     66                              );
     67
     68
     69/**
     70 * Plot the two luminosity functions created within psastroRefStarSubset
     71 * The luminosity functions are used to select a subset of reference stars,
     72 * so we plot the cutoff that defines this subset
     73 */
     74bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag,   ///< Log(n) for each magnitude bin
     75                                          psVector *Mag,     ///< magnitude bins
     76                                          pmLumFunc *lumFunc,///< Fit to the reference star luminosity function
     77                                          pmLumFunc *rawFunc ///< Fit to the raw star luminoisty function
     78                                           );
     79
     80
     81/**
     82 * Plot raw stars as determined from first pass astrometry fit
     83 * Called within psastroAstromGeuss
     84 */
     85bool pmAstromVisualPlotRawStars (psArray *rawstars, ///< Stars detected in the fpa
     86                                 pmFPA *fpa,  ///< structure describing the focal plane array
     87                                 pmChip *chip,  ///< structure describing the chip
     88                                 psMetadata *recipe ///< the recipe used in psastro
     89                                 );
     90
     91
     92/**
     93 * plot the location of references stars over the entire fpa
     94 * invoked during psastroChooseRefStars
     95 */
     96bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe);
     97
     98
     99/**
     100 * Plot the stars in a region, and indicate which stars are part of 'clumps'
     101 * These stars are flagged during astrometric fitting, since dense regions are
     102 * harder to cross-match than sparse ones. Called during psastroRemoveClumps.
     103 */
     104bool pmAstromVisualPlotRemoveClumps (psArray *input, ///< Array containing the field stars
     105                                    psImage *count, ///< A 2D histogram of the field star distribution
     106                                    int scale,      ///< The pixel size of the histogram
     107                                    float limit     ///< The minimum numuber of stars in a bin flagged as a clump
     108                                     );
     109
     110/**
     111 * Plots the chip corners in the FP before and after chips with inconsistent solutions have been fixed.
     112 * Invoked during psastroFixChips
     113 */
     114bool pmAstromVisualPlotFixChips (pmFPAfile *input, ///< focal plane array file
     115                                 psVector *xOld, ///< old X location of chip cornerss
     116                                 psVector *yOld ///< old Y location of chip corners
     117                                 );
     118
     119
     120/**
     121 * Assess the goodness of fit for a signle chip by
     122 * plotting the fit residuals
     123 * invoked during psastroOneChipFit
     124 */
     125bool pmAstromVisualPlotOneChipFit (psArray *rawstars, ///< stars detected in the image
     126                                  psArray *refstars, ///< reference stars over the same region
     127                                  psArray *match,    ///< contains which rawstars match to which refstars
     128                                  psMetadata *recipe ///< data reduction recipe
     129                                   );
     130
     131/**
     132 *  Plots the fpa chip corners projected on to the tangential plane before and after
     133 *  the astrometry solution has been applied. In psastroAstromGuessCheck, the old corners
     134 *  are then fit to the new corners to get a sense at how far off the initial WCS info was
     135 *  in offset, rotation, and scale. This procedure also plots the residuals of the fit from
     136 *  old to new coordinates
     137 */
     138bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, ///< P coordinates of chip corners before fitting
     139                                        psVector *cornerQo, ///< Q coordinates of chip corners before fitting
     140                                        psVector *cornerPn, ///< P coordinates of chip corners after fitting
     141                                        psVector *cornerQn, ///< Q coordinates of chip corners after fitting
     142                                        psVector *cornerPd, ///< P coordinate residuals of fit from old to new coordinates
     143                                        psVector *cornerQd  ///< Q coordinate residuals of fit from old to new coordinates
     144                                         );
     145
     146
     147/**
     148 *   plot the residuals between raw stars and ref stars after
     149 *   fitting in psastroMosaicOneChip
     150 */
     151bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) ;
     152
     153
     154/**
     155 * Plots the pixel scales of the fpa before they are
     156 * equalized in psastroMosaicCommonScale
     157 */
     158bool pmAstromVisualPlotCommonScale (pmFPA *fpa,         ///< the fpa
     159                                   psVector *oldScale  ///< the old pixel scale of each chip in the fpa
     160                                    );
     161
     162
     163/** pmAstromVisualPlotMosaicMatches
     164 * Plot the matches between raw and reference stars during pmAstromVisualMosaicSetMatch
     165 */
     166bool pmAstromVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe);
     167
     168
     169#endif

  • branches/cnb_branch_20090113/psModules/src/extras/Makefile.am

    r10610 r21208  
    77        psPipe.c \
    88        psIOBuffer.c \
    9         pmKapaPlots.c
     9        pmKapaPlots.c \
     10        pmVisual.c
    1011
    1112pkginclude_HEADERS = \
    … …  
    1314        psPipe.h \
    1415        psIOBuffer.h \
    15         pmKapaPlots.h
     16        pmKapaPlots.h \
     17        pmVisual.h
    1618
    1719CLEANFILES = *~

  • branches/cnb_branch_20090113/psModules/src/imcombine/Makefile.am

    r20049 r21208  
    1818        pmSubtractionStamps.c   \
    1919        pmSubtractionThreads.c  \
    20         pmPSFEnvelope.c
     20        pmPSFEnvelope.c         \
     21        pmSubtractionVisual.c
    2122
    2223pkginclude_HEADERS = \
    … …  
    3536        pmSubtractionStamps.h   \
    3637        pmSubtractionThreads.h  \
    37         pmPSFEnvelope.h
     38        pmPSFEnvelope.h         \
     39        pmSubtractionVisual.h
    3840
    3941CLEANFILES = *~

  • branches/cnb_branch_20090113/psModules/src/imcombine/pmSubtractionAnalysis.c

    r20052 r21208  
    1212
    1313#include "pmSubtractionAnalysis.h"
     14#include "pmSubtractionVisual.h"
    1415
    1516#define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
    … …  
    101102        }
    102103
     104        pmSubtractionVisualPlotConvKernels(convKernels);
    103105        psMetadataAddImage(analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE",
    104106                           PS_META_DUPLICATE_OK, "Realisations of kernel", convKernels);

  • branches/cnb_branch_20090113/psModules/src/imcombine/pmSubtractionEquation.c

    r20561 r21208  
    1313
    1414#include "pmSubtractionEquation.h"
    15 
     15#include "pmSubtractionVisual.h"
    1616
    1717//#define TESTING
    … …  
    695695    }
    696696
     697    pmSubtractionVisualPlotLeastSquares(stamps);
     698
    697699    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
    698700             psTimerClear("pmSubtractionCalculateEquation"));
    … …  
    10011003     }
    10021004
     1005    pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
    10031006    return true;
    10041007}

  • branches/cnb_branch_20090113/psModules/src/imcombine/pmSubtractionMatch.c

    r20718 r21208  
    1919#include "pmSubtractionThreads.h"
    2020#include "pmSubtractionMatch.h"
    21 
     21#include "pmSubtractionVisual.h"
    2222
    2323#define BG_STAT PS_STAT_ROBUST_MEDIAN   // Statistic to use for background
    … …  
    8686
    8787    memCheck("   extract stamps");
    88 
     88    pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1);
    8989    return true;
    9090}

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

    r20949 r21208  
    99#include <psVectorBracket.h>
    1010#include <pmKapaPlots.h>
     11#include <pmVisual.h>
    1112
    1213// XXX the following headers define constructs needed by the elements below
    … …  
    9899#include <pmSubtractionEquation.h>
    99100#include <pmReadoutCombine.h>
     101#include <pmSubtractionVisual.h>
    100102
    101103// the following headers are from psModule:objects

  • branches/cnb_branch_20090113/psastro/src/Makefile.am

    r20805 r21208  
    5353        psastroErrorCodes.c         \
    5454        psastroVersion.c            \
    55         psastroVisual.c             \
    5655        psastroDefineFiles.c        \
    5756        psastroAnalysis.c           \

  • branches/cnb_branch_20090113/psastro/src/psastro.h

    r20805 r21208  
    1414# define SIGN(X)  (((X) == 0) ? 0 : ((fabs((double)(X))) / (X)))
    1515
     16#if 0
    1617// this structure represents a fit to the logN / logS curve for a set of stars
    1718// logN = offset + slope * logS
    … …  
    2526    int sPeak;                          // sum of stars to peak bin
    2627} pmLumFunc;
     28#endif
    2729
    2830bool              psastroDataSave (pmConfig *config);
    … …  
    7779psString          psastroVersionLong(void);
    7880
    79 // psastroVisual functions
    80 bool psastroSetVisual (bool mode);
    81 bool psastroVisualClose();
    82 bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc);
    83 bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);
    84 bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe);
    85 bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit);
    86 bool psastroVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld);
    87 bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
    88 bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd);
    89 bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
    90 bool psastroVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale);
    91 bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe);
    92 
    9381// demo plots
    9482bool              psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);

  • branches/cnb_branch_20090113/psastro/src/psastroArguments.c

    r20805 r21208  
    8888    if ((N = psArgumentGet (argc, argv, "-visual"))) {
    8989        psArgumentRemove (N, &argc, argv);
    90         psastroSetVisual (true);
    9190        pmAstromSetVisual (true);
    9291    }

  • branches/cnb_branch_20090113/psastro/src/psastroAstromGuess.c

    r20805 r21208  
    2929    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    3030    if (!recipe) {
    31         psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
    32         return false;
     31        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
     32        return false;
    3333    }
    3434
    … …  
    3636    bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL");
    3737    if (!status) {
    38         useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
     38        useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
    3939    }
    4040
    … …  
    4242    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    4343    if (!input) {
    44         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    45         return false;
     44        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     45        return false;
    4646    }
    4747
    … …  
    4949    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    5050    if (!status) {
    51         psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
    52         return false;
    53     } 
     51        psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
     52        return false;
     53    }
    5454
    5555    psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32);
    … …  
    6767    bool bilevelAstrometry = false;
    6868    if (!useModel) {
    69         psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
     69        psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
    7070    }
    7171
    … …  
    7575        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    7676
    77         if (!useModel) {
    78             if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
    79         }
     77        if (!useModel) {
     78            if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
     79        }
    8080
    8181        if (newFPA) {
    8282            newFPA = false;
    83             while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
    84             while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
     83            while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
     84            while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
    8585            RAminSky = fpa->toSky->R - M_PI;
    8686            RAmaxSky = fpa->toSky->R + M_PI;
    8787        }
    8888
    89         // report and save the current best guess for the chip 0,0 pixel coordinates
    90         {
    91             psPlane ptCH, ptFP, ptTP;
    92             psSphere ptSky;
    93 
    94             ptCH.x = 0;
    95             ptCH.y = 0;
    96             psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    97             psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    98             psDeproject (&ptSky, &ptTP, fpa->toSky);
    99             psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    100 
    101             psVectorAppend (cornerL, ptFP.x);
    102             psVectorAppend (cornerM, ptFP.y);
    103             psVectorAppend (cornerP, ptTP.x);
    104             psVectorAppend (cornerQ, ptTP.y);
    105             psVectorAppend (cornerR, ptSky.r);
    106             psVectorAppend (cornerD, ptSky.d);
    107         }
     89        // report and save the current best guess for the chip 0,0 pixel coordinates
     90        {
     91            psPlane ptCH, ptFP, ptTP;
     92            psSphere ptSky;
     93
     94            ptCH.x = 0;
     95            ptCH.y = 0;
     96            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     97            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     98            psDeproject (&ptSky, &ptTP, fpa->toSky);
     99            psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     100
     101            psVectorAppend (cornerL, ptFP.x);
     102            psVectorAppend (cornerM, ptFP.y);
     103            psVectorAppend (cornerP, ptTP.x);
     104            psVectorAppend (cornerQ, ptTP.y);
     105            psVectorAppend (cornerR, ptSky.r);
     106            psVectorAppend (cornerD, ptSky.d);
     107        }
    108108
    109109        // apply the new WCS guess data to all of the data in the readouts
    … …  
    119119                if (rawstars == NULL) { continue; }
    120120
    121                 *nStars += rawstars->n;
     121                *nStars += rawstars->n;
    122122                for (int i = 0; i < rawstars->n; i++) {
    123123                    pmAstromObj *raw = rawstars->data[i];
    … …  
    138138                }
    139139
    140                 // dump or plot the resulting projected positions
    141                 if (psTraceGetLevel("psastro.dump") > 0) {
    142                     psastroDumpRawstars (rawstars, fpa, chip);
    143                 }
    144 
    145                 psastroVisualPlotRawStars(rawstars, fpa, chip, recipe);
    146 
    147                 if (psTraceGetLevel("psastro.plot") > 0) {
    148                     psastroPlotRawstars (rawstars, fpa, chip, recipe);
    149                 }
     140                // dump or plot the resulting projected positions
     141                if (psTraceGetLevel("psastro.dump") > 0) {
     142                    psastroDumpRawstars (rawstars, fpa, chip);
     143                }
     144
     145                pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe);
     146
     147                if (psTraceGetLevel("psastro.plot") > 0) {
     148                    psastroPlotRawstars (rawstars, fpa, chip, recipe);
     149                }
    150150            }
    151151        }
    … …  
    157157    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR",  PS_META_REPLACE, "", *nStars);
    158158    if (*nStars == 0) {
    159         psLogMsg ("psastro", 2, "no sources available for astrometry\n");
    160         psFree (view);
    161         return true;
    162     }
    163 
    164     psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 
    165               DEG_RAD*RAmin, DEG_RAD*DECmin,
    166               DEG_RAD*RAmax, DEG_RAD*DECmax);
     159        psLogMsg ("psastro", 2, "no sources available for astrometry\n");
     160        psFree (view);
     161        return true;
     162    }
     163
     164    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n",
     165              DEG_RAD*RAmin, DEG_RAD*DECmin,
     166              DEG_RAD*RAmax, DEG_RAD*DECmax);
    167167
    168168    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", RAmin);
    … …  
    203203    pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    204204    if (bilevelAstrometry) {
    205         if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    206             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    207             return false;
    208         }
     205        if (!pmAstromReadBilevelChip (chip, hdu->header)) {
     206            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     207            return false;
     208        }
    209209    } else {
    210         if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
    211             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    212             return false;
    213         }
     210        if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
     211            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     212            return false;
     213        }
    214214    }
    215215    return true;
    … …  
    225225    // load mosaic-level astrometry?
    226226    if (phu) {
    227         char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
    228         if (ctype) {
    229             *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    230         }
     227        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
     228        if (ctype) {
     229            *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     230        }
    231231    }
    232232    if (*bilevelAstrometry) {
    233         pmAstromReadBilevelMosaic (fpa, phu->header);
    234     } 
     233        pmAstromReadBilevelMosaic (fpa, phu->header);
     234    }
    235235    psFree (view);
    236236    return true;
    … …  
    245245    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    246246    if (!input) {
    247         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    248         return false;
     247        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     248        return false;
    249249    }
    250250
    … …  
    277277        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    278278
    279         // XXX we are currently inconsistent with marking the good vs the bad data
    280         // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
    281         // probably wrong since it implies there is no data!
    282 
    283         // skip chips for which the astrometry failed (NASTRO == 0)
    284         if (!chip->cells->n) goto skip_chip;
    285         pmCell *cell = chip->cells->data[0];
    286         if (!cell) goto skip_chip;
    287 
    288         if (!cell->readouts->n) goto skip_chip;
    289         pmReadout *readout = cell->readouts->data[0];
    290         if (!readout) goto skip_chip;
    291 
    292         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    293         if (!updates) goto skip_chip;
    294        
    295         int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
    296         if (!nAstro) goto skip_chip;
    297 
    298         float astError = psMetadataLookupF32 (&status, updates, "CERROR");
    299         if (fabs(astError) < 1e-6) goto skip_chip;
    300 
    301         psPlane ptCH, ptFP, ptTP;
    302         psSphere ptSky;
    303 
    304         ptCH.x = 0;
    305         ptCH.y = 0;
    306         psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    307         psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    308         psDeproject (&ptSky, &ptTP, fpa->toSky);
    309         psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    310 
    311         // new corner locations based on the calibrated astrometry
    312         psVectorAppend (cornerLn, ptFP.x);
    313         psVectorAppend (cornerMn, ptFP.y);
    314         psVectorAppend (cornerPn, ptTP.x);
    315         psVectorAppend (cornerQn, ptTP.y);
    316         psVectorAppend (cornerRn, ptSky.r);
    317         psVectorAppend (cornerDn, ptSky.d);
    318         psVectorAppend (cornerMK, 0);
    319         continue;
     279        // XXX we are currently inconsistent with marking the good vs the bad data
     280        // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
     281        // probably wrong since it implies there is no data!
     282
     283        // skip chips for which the astrometry failed (NASTRO == 0)
     284        if (!chip->cells->n) goto skip_chip;
     285        pmCell *cell = chip->cells->data[0];
     286        if (!cell) goto skip_chip;
     287
     288        if (!cell->readouts->n) goto skip_chip;
     289        pmReadout *readout = cell->readouts->data[0];
     290        if (!readout) goto skip_chip;
     291
     292        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     293        if (!updates) goto skip_chip;
     294
     295        int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
     296        if (!nAstro) goto skip_chip;
     297
     298        float astError = psMetadataLookupF32 (&status, updates, "CERROR");
     299        if (fabs(astError) < 1e-6) goto skip_chip;
     300
     301        psPlane ptCH, ptFP, ptTP;
     302        psSphere ptSky;
     303
     304        ptCH.x = 0;
     305        ptCH.y = 0;
     306        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     307        psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     308        psDeproject (&ptSky, &ptTP, fpa->toSky);
     309        psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     310
     311        // new corner locations based on the calibrated astrometry
     312        psVectorAppend (cornerLn, ptFP.x);
     313        psVectorAppend (cornerMn, ptFP.y);
     314        psVectorAppend (cornerPn, ptTP.x);
     315        psVectorAppend (cornerQn, ptTP.y);
     316        psVectorAppend (cornerRn, ptSky.r);
     317        psVectorAppend (cornerDn, ptSky.d);
     318        psVectorAppend (cornerMK, 0);
     319        continue;
    320320
    321321    skip_chip:
    322         // new corner locations based on the calibrated astrometry
    323         psVectorAppend (cornerLn, 0.0);
    324         psVectorAppend (cornerMn, 0.0);
    325         psVectorAppend (cornerPn, 0.0);
    326         psVectorAppend (cornerQn, 0.0);
    327         psVectorAppend (cornerRn, 0.0);
    328         psVectorAppend (cornerDn, 0.0);
    329         psVectorAppend (cornerMK, 1);
     322        // new corner locations based on the calibrated astrometry
     323        psVectorAppend (cornerLn, 0.0);
     324        psVectorAppend (cornerMn, 0.0);
     325        psVectorAppend (cornerPn, 0.0);
     326        psVectorAppend (cornerQn, 0.0);
     327        psVectorAppend (cornerRn, 0.0);
     328        psVectorAppend (cornerDn, 0.0);
     329        psVectorAppend (cornerMK, 1);
    330330    }
    331331
    … …  
    336336
    337337    for (int i = 0; i < cornerRo->n; i++) {
    338        
    339         psPlane ptTP;
    340         psSphere ptSky;
    341 
    342         ptSky.r = cornerRo->data.F32[i];
    343         ptSky.d = cornerDo->data.F32[i];
    344 
    345         psProject (&ptTP, &ptSky, fpa->toSky);
    346         psVectorAppend (cornerPs, ptTP.x);
    347         psVectorAppend (cornerQs, ptTP.y);
     338
     339        psPlane ptTP;
     340        psSphere ptSky;
     341
     342        ptSky.r = cornerRo->data.F32[i];
     343        ptSky.d = cornerDo->data.F32[i];
     344
     345        psProject (&ptTP, &ptSky, fpa->toSky);
     346        psVectorAppend (cornerPs, ptTP.x);
     347        psVectorAppend (cornerQs, ptTP.y);
    348348    }
    349349
    … …  
    351351    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
    352352    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
    353    
     353
    354354    // fit the valid chips, mask the invalid chips
    355355    psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs);
    356356    psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs);
    357    
     357
    358358    // apply the linear fit...
    359359    psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs);
    … …  
    364364    psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf);
    365365
    366     psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);
     366    pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);
    367367
    368368    psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    … …  
    374374    float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);
    375375    float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]);
    376    
     376
    377377    psLogMsg ("psastro", 3, "boresite offset  : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]);
    378378    psLogMsg ("psastro", 3, "boresite angle   : %f, scale: %f", angle*PS_DEG_RAD, scale);
    … …  
    382382    psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER");
    383383    if (!header) {
    384         header = psMetadataAlloc();
    385         psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
    386         psFree (header);  // drop this reference
     384        header = psMetadataAlloc();
     385        psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
     386        psFree (header);  // drop this reference
    387387    }
    388388
    … …  
    395395
    396396    if (DEBUG) {
    397         FILE *f = fopen ("corners.dat", "w");
    398         for (int i = 0; i < cornerRo->n; i++) {
    399             fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
    400                      cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
    401                      cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
    402         }
    403         fclose (f);
     397        FILE *f = fopen ("corners.dat", "w");
     398        for (int i = 0; i < cornerRo->n; i++) {
     399            fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
     400                     cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
     401                     cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
     402        }
     403        fclose (f);
    404404    }
    405405
    … …  
    424424    psFree (map);
    425425    psFree (view);
    426    
     426
    427427
    428428    return true;

  • branches/cnb_branch_20090113/psastro/src/psastroCleanup.c

    r20805 r21208  
    44
    55    psFree (config);
     6    pmAstromVisualClose ();
    67
    78    psTimerStop ();
    … …  
    1112    pmConceptsDone ();
    1213    pmConfigDone ();
    13     psastroVisualClose ();
    14     pmAstromVisualClose ();
    1514    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
    1615    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");

  • branches/cnb_branch_20090113/psastro/src/psastroFixChips.c

    r20805 r21208  
    1313    bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
    1414    if (!status) {
    15         fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
     15        fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
    1616    }
    1717    if (!fixChips) return true;
    … …  
    2525    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    2626    if (!input) {
    27         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    28         return false;
     27        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     28        return false;
    2929    }
    3030
    … …  
    4848    // files associated with the science image
    4949    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    50         psError (PS_ERR_IO, false, "Can't load the astrometry model file");
    51         return false;
     50        psError (PS_ERR_IO, false, "Can't load the astrometry model file");
     51        return false;
    5252    }
    5353
    … …  
    6464
    6565    if (DEBUG) {
    66         f = fopen ("corners.raw.dat", "w");
    67         chipName = NULL;
     66        f = fopen ("corners.raw.dat", "w");
     67        chipName = NULL;
    6868    }
    6969
    7070    pmChip *obsChip = NULL;
    7171    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    72         if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    73 
    74         // XXX we are currently inconsistent with marking the good vs the bad data
    75         // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
    76         // probably wrong since it implies there is no data!
    77 
    78         // skip chips for which the astrometry failed (NASTRO == 0)
    79         if (!obsChip->cells->n) continue;
    80         pmCell *cell = obsChip->cells->data[0];
    81         if (!cell) continue;
    82 
    83         if (!cell->readouts->n) continue;
    84         pmReadout *readout = cell->readouts->data[0];
    85         if (!readout) continue;
    86 
    87         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    88         if (!updates) continue;
    89        
    90         int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
    91         if (!nAstro) continue;
    92 
    93         // set the chip astrometry using the astrom file
    94         pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    95 
    96         psRegion *region = pmChipPixels (obsChip);
    97         psPlane ptCP, ptFP;
    98 
    99         ptCP.x = region->x0; ptCP.y = region->y0;
    100         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    101         xObs->data.F32[nPts] = ptFP.x;
    102         yObs->data.F32[nPts] = ptFP.y;
    103         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    104         xRef->data.F32[nPts] = ptFP.x;
    105         yRef->data.F32[nPts] = ptFP.y;
    106 
    107         if (DEBUG) {
    108             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    109             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    110         }
    111         nPts ++;
    112 
    113         ptCP.x = region->x0; ptCP.y = region->y1;
    114         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    115         xObs->data.F32[nPts] = ptFP.x;
    116         yObs->data.F32[nPts] = ptFP.y;
    117         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    118         xRef->data.F32[nPts] = ptFP.x;
    119         yRef->data.F32[nPts] = ptFP.y;
    120 
    121         if (DEBUG) {
    122             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    123             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    124         }
    125         nPts ++;
    126 
    127         ptCP.x = region->x1; ptCP.y = region->y1;
    128         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    129         xObs->data.F32[nPts] = ptFP.x;
    130         yObs->data.F32[nPts] = ptFP.y;
    131         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    132         xRef->data.F32[nPts] = ptFP.x;
    133         yRef->data.F32[nPts] = ptFP.y;
    134 
    135         if (DEBUG) {
    136             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    137             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    138         }
    139         nPts ++;
    140 
    141         ptCP.x = region->x1; ptCP.y = region->y0;
    142         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    143         xObs->data.F32[nPts] = ptFP.x;
    144         yObs->data.F32[nPts] = ptFP.y;
    145         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    146         xRef->data.F32[nPts] = ptFP.x;
    147         yRef->data.F32[nPts] = ptFP.y;
    148 
    149         if (DEBUG) {
    150             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    151             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    152         }
    153         nPts ++;
    154 
    155         psFree (region);
     72        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     73
     74        // XXX we are currently inconsistent with marking the good vs the bad data
     75        // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
     76        // probably wrong since it implies there is no data!
     77
     78        // skip chips for which the astrometry failed (NASTRO == 0)
     79        if (!obsChip->cells->n) continue;
     80        pmCell *cell = obsChip->cells->data[0];
     81        if (!cell) continue;
     82
     83        if (!cell->readouts->n) continue;
     84        pmReadout *readout = cell->readouts->data[0];
     85        if (!readout) continue;
     86
     87        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     88        if (!updates) continue;
     89
     90        int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
     91        if (!nAstro) continue;
     92
     93        // set the chip astrometry using the astrom file
     94        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     95
     96        psRegion *region = pmChipPixels (obsChip);
     97        psPlane ptCP, ptFP;
     98
     99        ptCP.x = region->x0; ptCP.y = region->y0;
     100        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     101        xObs->data.F32[nPts] = ptFP.x;
     102        yObs->data.F32[nPts] = ptFP.y;
     103        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     104        xRef->data.F32[nPts] = ptFP.x;
     105        yRef->data.F32[nPts] = ptFP.y;
     106
     107        if (DEBUG) {
     108            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     109            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     110        }
     111        nPts ++;
     112
     113        ptCP.x = region->x0; ptCP.y = region->y1;
     114        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     115        xObs->data.F32[nPts] = ptFP.x;
     116        yObs->data.F32[nPts] = ptFP.y;
     117        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     118        xRef->data.F32[nPts] = ptFP.x;
     119        yRef->data.F32[nPts] = ptFP.y;
     120
     121        if (DEBUG) {
     122            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     123            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     124        }
     125        nPts ++;
     126
     127        ptCP.x = region->x1; ptCP.y = region->y1;
     128        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     129        xObs->data.F32[nPts] = ptFP.x;
     130        yObs->data.F32[nPts] = ptFP.y;
     131        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     132        xRef->data.F32[nPts] = ptFP.x;
     133        yRef->data.F32[nPts] = ptFP.y;
     134
     135        if (DEBUG) {
     136            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     137            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     138        }
     139        nPts ++;
     140
     141        ptCP.x = region->x1; ptCP.y = region->y0;
     142        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     143        xObs->data.F32[nPts] = ptFP.x;
     144        yObs->data.F32[nPts] = ptFP.y;
     145        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     146        xRef->data.F32[nPts] = ptFP.x;
     147        yRef->data.F32[nPts] = ptFP.y;
     148
     149        if (DEBUG) {
     150            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     151            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     152        }
     153        nPts ++;
     154
     155        psFree (region);
    156156    }
    157157    xObs->n = yObs->n = xRef->n = yRef->n = nPts;
    158158    if (DEBUG) fclose (f);
    159        
     159
    160160    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
    161  
     161
    162162    psVector *mask = psVectorAlloc (nPts, PS_TYPE_U8);
    163163    psVectorInit (mask, 0);
    … …  
    167167
    168168    for (int i = 0; i < 3; i++) {
    169         psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
    170         psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
    171 
    172         psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
    173         psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
     169        psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
     170        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
     171
     172        psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
     173        psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
    174174    }
    175175
    176176    // loop over all chips, select the outliers, and replace the measured astrometry with the model
    177     // the measured transformation above must be applied to make the comparison, and also then applied to the 
     177    // the measured transformation above must be applied to make the comparison, and also then applied to the
    178178    // model transformation
    179179
    180180    if (DEBUG) {
    181         f = fopen ("corners.fit.dat", "w");
    182         for (int i = 0; i < xObs->n; i++) {
    183             psPlane obsCoord, refCoord;
    184             refCoord.x = xRef->data.F32[i];
    185             refCoord.y = yRef->data.F32[i];
    186             psPlaneTransformApply (&obsCoord, map, &refCoord);
    187             fprintf (f, "%f %f  %f %f  %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y);
    188         }
    189         fclose (f);
     181        f = fopen ("corners.fit.dat", "w");
     182        for (int i = 0; i < xObs->n; i++) {
     183            psPlane obsCoord, refCoord;
     184            refCoord.x = xRef->data.F32[i];
     185            refCoord.y = yRef->data.F32[i];
     186            psPlaneTransformApply (&obsCoord, map, &refCoord);
     187            fprintf (f, "%f %f  %f %f  %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y);
     188        }
     189        fclose (f);
    190190    }
    191191
    … …  
    199199
    200200    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    201         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
    202         if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    203 
    204         // set the chip astrometry using the astrom file
    205         pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    206 
    207         // bad Astrometry test:  ref pixel or angle outside nominal
    208 
    209         psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
    210         psPlane obsCoord, refCoord, tmpCoord;
    211 
    212         // find location of 0,0 pixel in focal plane coords for this chip
    213         psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
    214 
    215         // find location of 0,0 pixel in focal plane coords for ref chip
    216         // apply the global field rotation and offset before comparing
    217         psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
    218         psPlaneTransformApply (&refCoord, map, &tmpCoord);
    219    
    220         psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
    221         psPlane obsOffPt, refOffPt;
    222 
    223         // find location of 0,0 pixel in focal plane coords for this chip
    224         psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
    225 
    226         // find location of 0,0 pixel in focal plane coords for ref chip
    227         psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
    228         psPlaneTransformApply (&refOffPt, map, &tmpCoord);
    229    
    230         double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
    231         double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
    232 
    233         bool badAstrom = false;
    234         badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
    235         badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
    236         badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
    237 
    238         fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    239 
    240         // XXX for now, just use first readout
    241         pmCell *cell = obsChip->cells->data[0];
    242         pmReadout *readout = cell->readouts->data[0];
    243 
    244         // update the header (pull or create local view to entry on readout->analysis)
    245         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    246         if (!updates) {
    247             updates = psMetadataAlloc ();
    248             psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
    249             psFree (updates);
    250         }
    251 
    252         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x);
    253         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
    254         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
    255 
    256         // for successful chips, save the measured offsets in the header
    257         if (!badAstrom) continue;
    258 
    259         // XXX for now, let's just fail on the bad chips.  In the future, let's try to recover, but we still need to
    260         // catch the failures relative to the model
    261         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0);
    262         continue;
    263 
    264         psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
    265                   view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    266 
    267         psFree (obsChip->toFPA);
    268         psFree (obsChip->fromFPA);
    269 
    270         // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
    271         // XXX this only works if toTPA is at most a linear transformation
    272         psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
    273         for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
    274             for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
    275                 double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    276                 double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    277                 toFPA->x->coeff[i][j] = f1 + f2;
    278 
    279                 double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    280                 double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    281                 toFPA->y->coeff[i][j] = g1 + g2;
    282             }
    283         }
    284         toFPA->x->coeff[0][0] += map->x->coeff[0][0];
    285         toFPA->y->coeff[0][0] += map->y->coeff[0][0];
    286 
    287         psRegion *region = pmChipPixels (obsChip);
    288         obsChip->toFPA   = toFPA;
    289         obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
    290         psFree (region);
    291    
    292         // use the new position to re-try the match fit
    293         // select the raw objects for this readout
    294         psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    295         if (rawstars == NULL) { continue; }
    296 
    297         // select the raw objects for this readout
    298         psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    299         if (refstars == NULL) { continue; }
    300 
    301         // the absolute minimum number of stars is 4 (for order = 1)
    302         if ((rawstars->n < 4) || (refstars->n < 4)) {
    303             readout->data_exists = false;
    304             psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
    305                       rawstars->n, refstars->n);
    306             continue;
    307         }
    308 
    309         psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
    310 
    311         // XXX update the header with info to reflect the failure
    312         if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
    313             readout->data_exists = false;
    314             psLogMsg ("psastro", 3, "failed to find a solution\n");
    315             continue;
    316         }
    317 
    318         pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
    319     }
    320 
    321     psastroVisualPlotFixChips (input, xObs, yObs);
     201        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
     202        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     203
     204        // set the chip astrometry using the astrom file
     205        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     206
     207        // bad Astrometry test:  ref pixel or angle outside nominal
     208
     209        psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
     210        psPlane obsCoord, refCoord, tmpCoord;
     211
     212        // find location of 0,0 pixel in focal plane coords for this chip
     213        psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
     214
     215        // find location of 0,0 pixel in focal plane coords for ref chip
     216        // apply the global field rotation and offset before comparing
     217        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
     218        psPlaneTransformApply (&refCoord, map, &tmpCoord);
     219
     220        psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
     221        psPlane obsOffPt, refOffPt;
     222
     223        // find location of 0,0 pixel in focal plane coords for this chip
     224        psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
     225
     226        // find location of 0,0 pixel in focal plane coords for ref chip
     227        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
     228        psPlaneTransformApply (&refOffPt, map, &tmpCoord);
     229
     230        double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
     231        double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
     232
     233        bool badAstrom = false;
     234        badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
     235        badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
     236        badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
     237
     238        fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     239
     240        // XXX for now, just use first readout
     241        pmCell *cell = obsChip->cells->data[0];
     242        pmReadout *readout = cell->readouts->data[0];
     243
     244        // update the header (pull or create local view to entry on readout->analysis)
     245        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     246        if (!updates) {
     247            updates = psMetadataAlloc ();
     248            psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     249            psFree (updates);
     250        }
     251
     252        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x);
     253        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
     254        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
     255
     256        // for successful chips, save the measured offsets in the header
     257        if (!badAstrom) continue;
     258
     259        // XXX for now, let's just fail on the bad chips.  In the future, let's try to recover, but we still need to
     260        // catch the failures relative to the model
     261        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0);
     262        continue;
     263
     264        psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
     265                  view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     266
     267        psFree (obsChip->toFPA);
     268        psFree (obsChip->fromFPA);
     269
     270        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
     271        // XXX this only works if toTPA is at most a linear transformation
     272        psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
     273        for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
     274            for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
     275                double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     276                double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     277                toFPA->x->coeff[i][j] = f1 + f2;
     278
     279                double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     280                double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     281                toFPA->y->coeff[i][j] = g1 + g2;
     282            }
     283        }
     284        toFPA->x->coeff[0][0] += map->x->coeff[0][0];
     285        toFPA->y->coeff[0][0] += map->y->coeff[0][0];
     286
     287        psRegion *region = pmChipPixels (obsChip);
     288        obsChip->toFPA   = toFPA;
     289        obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
     290        psFree (region);
     291
     292        // use the new position to re-try the match fit
     293        // select the raw objects for this readout
     294        psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     295        if (rawstars == NULL) { continue; }
     296
     297        // select the raw objects for this readout
     298        psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     299        if (refstars == NULL) { continue; }
     300
     301        // the absolute minimum number of stars is 4 (for order = 1)
     302        if ((rawstars->n < 4) || (refstars->n < 4)) {
     303            readout->data_exists = false;
     304            psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
     305                      rawstars->n, refstars->n);
     306            continue;
     307        }
     308
     309        psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
     310
     311        // XXX update the header with info to reflect the failure
     312        if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
     313            readout->data_exists = false;
     314            psLogMsg ("psastro", 3, "failed to find a solution\n");
     315            continue;
     316        }
     317
     318        pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
     319    }
     320
     321    pmAstromVisualPlotFixChips (input, xObs, yObs);
    322322    psFree (xObs);
    323323    psFree (yObs);

  • branches/cnb_branch_20090113/psastro/src/psastroLoadRefstars.c

    r20805 r21208  
    138138    }
    139139
    140     psastroVisualPlotRefStars (refstars, recipe);
     140    pmAstromVisualPlotRefStars (refstars, recipe);
    141141
    142142    if (psTraceGetLevel("psastro.plot") > 0) {
    … …  
    239239    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    240240    if (!input) {
    241         psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data");
    242         photcode = psStringCopy ("NONE");
    243         return photcode;
     241        psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data");
     242        photcode = psStringCopy ("NONE");
     243        return photcode;
    244244    }
    245245    assert (input->fpa);
    … …  
    247247    *maxRho = psMetadataLookupF32(&status, recipe, "DVO.GETSTAR.MAX.RHO");
    248248    if (!status) {
    249         psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe");
    250         return NULL;
     249        psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe");
     250        return NULL;
    251251    }
    252252
    … …  
    266266    if (!status) ESCAPE ("missing DVO.GETSTAR.MIN.MAG.INST");
    267267
    268     // PHOTCODE.DATA is a multi of metadata items 
     268    // PHOTCODE.DATA is a multi of metadata items
    269269    psListIterator *iter = psListIteratorAlloc(item->data.list, PS_LIST_HEAD, false);
    270270
    271271    psMetadataItem *refItem = NULL;
    272272    while ((refItem = psListGetAndIncrement (iter))) {
    273         if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder");
    274    
    275         char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
    276         if (!status) {
    277             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    278             continue;
    279         }
    280 
    281         // does this entry match the current filter?
    282         if (strcmp (refFilter, filter)) continue;
    283 
    284         psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter);
    285 
    286         float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
    287         if (!status) {
    288             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    289             continue;
    290         }
    291         photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE");
    292         if (!status) {
    293             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    294             continue;
    295         }
    296 
    297         // convert the minInst to a calibrated minimum magnitude
    298         *minMag = minInst + 2.5*log10(exptime) + zeropt;
    299 
    300         psFree (iter);
    301         return photcode;
    302     }   
     273        if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder");
     274
     275        char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
     276        if (!status) {
     277            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     278            continue;
     279        }
     280
     281        // does this entry match the current filter?
     282        if (strcmp (refFilter, filter)) continue;
     283
     284        psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter);
     285
     286        float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
     287        if (!status) {
     288            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     289            continue;
     290        }
     291        photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE");
     292        if (!status) {
     293            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     294            continue;
     295        }
     296
     297        // convert the minInst to a calibrated minimum magnitude
     298        *minMag = minInst + 2.5*log10(exptime) + zeropt;
     299
     300        psFree (iter);
     301        return photcode;
     302    }
    303303    psFree (iter);
    304304
    … …  
    306306    photcode = psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR.PHOTCODE");
    307307    PS_ASSERT (photcode, NULL);
    308        
     308
    309309    // give up and use fixed value
    310310    *minMag = psMetadataLookupF32(NULL, recipe, "DVO.GETSTAR.MIN.MAG");

  • branches/cnb_branch_20090113/psastro/src/psastroLuminosityFunction.c

    r20805 r21208  
    129129    lumFunc->sPeak = sPeak;
    130130
    131     psastroVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc);
     131    pmAstromVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc);
    132132
    133133    psFree (lnMag);

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

    r20803 r21208  
    148148
    149149    //plot results
    150     psastroVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);
     150    pmAstromVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);
    151151
    152152    psFree (fitStats);

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

    r20805 r21208  
    5757                psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
    5858
    59                 psastroVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);
     59                pmAstromVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);
    6060
    6161                // XXX drop the old one

  • branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c

    r20805 r21208  
    1111
    1212    // default value for match/fit : radius is in pixels
    13     REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
     13    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32);
    1414
    1515    // run the match/fit sequence NITER times
    16     REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 
     16    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32);
    1717
    1818    // correct radius to FP units (physical pixel scale in microns per pixel)
    19     REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 
     19    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32);
    2020    RADIUS *= pixelScale;
    2121
    … …  
    3232
    3333    for (int iter = 0; iter < nIter; iter++) {
    34        
    35         char name[128];
    3634
    37         sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
    38         float radius = psMetadataLookupF32 (&status, recipe, name);
    39         radius *= pixelScale;
    40         if (!status || (radius == 0.0)) {
    41             radius = RADIUS;
    42         }
     35        char name[128];
     36
     37        sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
     38        float radius = psMetadataLookupF32 (&status, recipe, name);
     39        radius *= pixelScale;
     40        if (!status || (radius == 0.0)) {
     41            radius = RADIUS;
     42        }
    4343
    4444
    45         // use small radius to match stars
    46         match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
    47         if (match == NULL) {
    48             psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
    49             return false;
    50         }
     45        // use small radius to match stars
     46        match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
     47        if (match == NULL) {
     48            psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
     49            return false;
     50        }
    5151
    52         // modify the order to correspond to the actual number of matched stars:
    53         int Ndof_min = 3;
    54         int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
    55         order = PS_MIN (order, order_max);
     52        // modify the order to correspond to the actual number of matched stars:
     53        int Ndof_min = 3;
     54        int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
     55        order = PS_MIN (order, order_max);
    5656
    57         // if ((match->n < 11) && (order >= 3)) order = 2;
    58         // if ((match->n <  7) && (order >= 2)) order = 1;
    59         // if ((match->n <  4) && (order >= 1)) order = 0;
     57        // if ((match->n < 11) && (order >= 3)) order = 2;
     58        // if ((match->n <  7) && (order >= 2)) order = 1;
     59        // if ((match->n <  4) && (order >= 1)) order = 0;
    6060
    61         if (order < 1) {
    62             psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
    63             psFree (match);
    64             return false;
    65         }
     61        if (order < 1) {
     62            psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
     63            psFree (match);
     64            return false;
     65        }
    6666
    67         // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
    68         psFree (chip->toFPA);
    69         chip->toFPA = psPlaneTransformAlloc (order, order);
    70         for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    71             for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    72                 if (i + j > order) {
    73                     chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
    74                     chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
    75                 }
    76             }
    77         }
     67        // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
     68        psFree (chip->toFPA);
     69        chip->toFPA = psPlaneTransformAlloc (order, order);
     70        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     71            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     72                if (i + j > order) {
     73                    chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
     74                    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
     75                }
     76            }
     77        }
    7878
    79         // XXX allow statitic to be set by the user
    80         // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    81         fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    82         fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
    83         fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
     79        // XXX allow statitic to be set by the user
     80        // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     81        fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     82        fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
     83        fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
    8484
    85         // improved fit for astrometric terms
    86         results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
    87         if (!results) {
    88             psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
    89             psFree (match);
    90             psFree (fitStats);
    91             return false;
    92         }
    93    
    94         // determine fromFPA transformation and apply new transformation to raw & ref stars
    95         psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    96    
    97         // toSky converts from FPA & TPA units (microns) to sky units (radians)
    98         float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    99         // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
    100         float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
    101         int astNstar = results->yStats->clippedNvalues;
    102         psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     85        // improved fit for astrometric terms
     86        results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
     87        if (!results) {
     88            psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
     89            psFree (match);
     90            psFree (fitStats);
     91            return false;
     92        }
    10393
    104         if (iter < nIter - 1) {
    105             psFree (fitStats);
    106             psFree (results);
    107             psFree (match);
    108         }
     94        // determine fromFPA transformation and apply new transformation to raw & ref stars
     95        psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     96
     97        // toSky converts from FPA & TPA units (microns) to sky units (radians)
     98        float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
     99        // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
     100        float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     101        int astNstar = results->yStats->clippedNvalues;
     102        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     103
     104        if (iter < nIter - 1) {
     105            psFree (fitStats);
     106            psFree (results);
     107            psFree (match);
     108        }
    109109    }
    110110
    … …  
    127127    if (astError > maxError) {
    128128        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    129         validSolution = false;
     129        validSolution = false;
    130130    }
    131131    if (astNstar < minNstar) {
    132132        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
    133         validSolution = false;
     133        validSolution = false;
    134134    }
    135135
    … …  
    138138    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
    139139    if (validSolution) {
    140         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
    141         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     140        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
     141        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
    142142    } else {
    143         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
    144         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     143        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
     144        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
    145145    }
    146146    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
    … …  
    148148    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
    149149    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    150    
     150
    151151    // XXX check if we correctly applied the new transformation:
    152152    if (psTraceGetLevel("psastro.dump") > 0) {
    153         psastroDumpRawstars (rawstars, fpa, chip);
    154         psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
    155         psastroDumpStars (refstars, "refstars.cal.dat");
     153        psastroDumpRawstars (rawstars, fpa, chip);
     154        psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
     155        psastroDumpStars (refstars, "refstars.cal.dat");
    156156    }
    157157
    158     psastroVisualPlotOneChipFit (rawstars, refstars, match, recipe);
     158    pmAstromVisualPlotOneChipFit (rawstars, refstars, match, recipe);
    159159
    160160    if (psTraceGetLevel("psastro.plot") > 0) {
    161         psastroPlotOneChipFit (rawstars, refstars, match, recipe);
     161        psastroPlotOneChipFit (rawstars, refstars, match, recipe);
    162162    }
    163163

  • branches/cnb_branch_20090113/psastro/src/psastroRemoveClumps.c

    r20805 r21208  
    5555    psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit);
    5656
    57     psastroVisualPlotRemoveClumps (input, count, scale, limit);
     57    pmAstromVisualPlotRemoveClumps (input, count, scale, limit);
    5858
    5959    // find and exclude objects in bad pixels

  • branches/cnb_branch_20090113/psastro/src/psastroUtils.c

    r20804 r21208  
    4242            if (!chip->toFPA) { continue; }
    4343
    44             if (chip->cells->n == 0) { continue; }
    45             pmCell *cell = chip->cells->data[0];
     44            if (chip->cells->n == 0) { continue; }
     45            pmCell *cell = chip->cells->data[0];
    4646            if (!cell->process || !cell->file_exists) { continue; }
    4747
    48             if (cell->readouts->n == 0) { continue; }
    49             pmReadout *readout = cell->readouts->data[0];
    50             if (! readout->data_exists) { continue; }
     48            if (cell->readouts->n == 0) { continue; }
     49            pmReadout *readout = cell->readouts->data[0];
     50            if (! readout->data_exists) { continue; }
    5151
    5252            pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]);
    … …  
    8888    psastroMosaicSetAstrom (fpa);
    8989    if (!useExternal) {
    90         psastroVisualPlotCommonScale (fpa, oldScale);
     90        pmAstromVisualPlotCommonScale (fpa, oldScale);
    9191    }
    9292    psFree (oldScale);
Note: See TracChangeset for help on using the changeset viewer.