IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15665


Ignore:
Timestamp:
Nov 20, 2007, 9:02:55 PM (18 years ago)
Author:
eugene
Message:

adding a re-scale option to pmAstrometryObjects grid match; working on astrometry I/O tables

Location:
trunk/psModules/src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/psModules/src/astrom/pmAstrometryObjects.c

    r15102 r15665  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2007-09-29 02:28:04 $
     10*  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2007-11-21 07:02:55 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    269269
    270270/******************************************************************************
    271 pmAstromRotateObj(old, center, angle): rotate the focal-plane coordinates
     271pmAstromRotateObj(old, center, angle, angle): rotate & scale the focal-plane coordinates
    272272about the center coordinate angle specified in radians
    273273 ******************************************************************************/
     
    275275    const psArray *old,
    276276    psPlane center,
    277     double angle)
     277    double angle,
     278    double scale)
    278279{
    279280    PS_ASSERT_PTR_NON_NULL(old, NULL);
     
    284285
    285286    psArray *new = psArrayAlloc (old->n);
    286     double cs = cos(angle);
    287     double sn = sin(angle);
     287    double cs = scale*cos(angle);
     288    double sn = scale*sin(angle);
    288289    double xCenter = center.x;
    289290    double yCenter = center.y;
     
    326327    //    stats->offset = {0, 0, 0, 0};
    327328    stats->angle     = 0.0;
     329    stats->scale     = 1.0;
    328330    stats->minMetric = 0.0;
    329331    stats->minVar    = 0.0;
     
    703705    center.y = 0.5*(yMin + yMax);
    704706
     707    double minScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.SCALE");
     708    double maxScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.SCALE");
     709    double delScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.SCALE");
     710
    705711    double minAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE");
    706712    double maxAngle = RAD_DEG*psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE");
     
    709715
    710716    minStat->minMetric = 1e10;
    711     for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
    712         rot = pmAstromRotateObj (raw, center, angle);
    713         newStat = pmAstromGridAngle (rot, ref, config);
    714         newStat->angle  = angle;
    715         newStat->center = center;
    716 
    717         if (newStat->minMetric < minStat->minMetric) {
    718             *minStat = *newStat;
    719             psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *",
    720                       minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
    721         } else {
    722             psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)",
    723                       newStat->offset.x, newStat->offset.y, DEG_RAD*newStat->angle, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric));
    724 
    725         }
    726         psFree (newStat);
    727         psFree (rot);
    728     }
    729     psLogMsg ("psModule.astrom.grid.match", 4, "grid best - offset: %7.2f,%7.2f @ %6.1f deg (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)",
    730               minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
     717    for (double scale = minScale; scale <= maxScale; scale += delScale) {
     718        for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
     719            rot = pmAstromRotateObj (raw, center, angle, scale);
     720            newStat = pmAstromGridAngle (rot, ref, config);
     721            newStat->angle  = angle;
     722            newStat->scale  = scale;
     723            newStat->center = center;
     724
     725            if (newStat->minMetric < minStat->minMetric) {
     726                *minStat = *newStat;
     727                psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *",
     728                          minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
     729            } else {
     730                psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)",
     731                          newStat->offset.x, newStat->offset.y, DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric));
     732
     733            }
     734            psFree (newStat);
     735            psFree (rot);
     736        }
     737    }
     738    psLogMsg ("psModule.astrom.grid.match", 4, "grid best - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)",
     739              minStat->offset.x, minStat->offset.y, DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
    731740
    732741    // I need to decide if a solution is likely to be a good solution or just a mis-match
     
    756765    int nBin, xBin, yBin;
    757766
    758     rot = pmAstromRotateObj (raw, stats->center, stats->angle);
     767    rot = pmAstromRotateObj (raw, stats->center, stats->angle, stats->scale);
    759768
    760769    // sampling scale of the grid
     
    846855    PS_ASSERT_POLY_NON_NULL(map->y, NULL);
    847856
    848     double cs = cos (stat->angle);
    849     double sn = sin (stat->angle);
     857    double cs = stat->scale * cos (stat->angle);
     858    double sn = stat->scale * sin (stat->angle);
    850859
    851860    double dx = (map->x->coeff[0][0] - stat->center.x);
  • trunk/psModules/src/astrom/pmAstrometryObjects.h

    r14649 r15665  
    44 * @author EAM, IfA
    55 *
    6  * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-08-23 23:44:05 $
     6 * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2007-11-21 07:02:55 $
    88 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    99 */
     
    6868    psPlane center;                     ///<
    6969    psPlane offset;                     ///<
     70    double  scale;                      ///<
    7071    double  angle;                      ///<
    7172    double  minMetric;                  ///<
     
    142143    const psArray *old,
    143144    psPlane center,
    144     double angle
     145    double angle,
     146    double scale
    145147);
    146148
  • trunk/psModules/src/astrom/pmAstrometryTable.c

    r15562 r15665  
    88 *  @author EAM, IfA
    99 *
    10  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2007-11-10 01:09:20 $
     10 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2007-11-21 07:02:55 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    3636#include "pmFPAview.h"
    3737#include "pmFPAfile.h"
     38#include "pmFPAExtent.h"
     39#include "pmAstrometryWCS.h"
     40#include "pmAstrometryRegions.h"
    3841#include "pmAstrometryTable.h"
    39 
    4042
    4143/********************* CheckDataStatus functions *****************************/
     
    181183
    182184// write out chip-level Astrometry data for this chip
     185// XXX this function writes out astrometry elements for a chip->fpa->tp->sky layer set
    183186// XXX we need to ensure the chip and fpa level data are written
    184187// XXX can we write the rows one at a time (or groups at a time?)
     188// XXX allow this to be chip level, or just force to be fpa?  is it valid at the chip level??
    185189bool pmAstromWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    186190{
    187     bool status;
    188     // *** define the EXTNAME values used for the table data, and residual image segments ***
    189 
    190     // lookup the EXTNAME values used for table data and image header segments
    191     char *rule = NULL;
    192 
    193     // Menu of EXTNAME rules
    194     psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
    195     if (!menu) {
    196         psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
    197         return false;
    198     }
    199 
    200     // EXTNAME for image header
    201     rule = psMetadataLookupStr(&status, menu, "ASTROMETRY");
    202     if (!rule) {
    203         psError(PS_ERR_UNKNOWN, false, "missing entry for ASTROMETRY in EXTNAME.RULES in camera.config");
    204         return false;
    205     }
    206 
    207     // XXX need to finish this: uncomment when used
    208     # if (0)
    209     char *extname = pmFPAfileNameFromRule (rule, file, view);
    210     # endif
    211 
    212     // write the chip elements in the following form:
    213     // chipID, ...?  I forget the form: look at ICD
    214 
     191    psArray *table;
     192    psMetadata *row;
     193    psMetadata *header;
     194
     195    // write out the blank PHU to start the file
     196    psFitsWriteBlank (file->fits, NULL, NULL);
     197
     198    // first layer is the sky
     199    header = psMetadataAlloc();
     200    psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "SKY");
     201    psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "NONE");
     202    psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "NONE");
     203    psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "NONE");
     204
     205    table = psArrayAllocEmpty (1);
     206    row = psMetadataAlloc ();
     207    psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "SKY");
     208    psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "NONE");
     209   
     210    psArrayAdd (table, 100, row);
     211    psFree (row);
     212
     213    if (!psFitsWriteTable (file->fits, header, table, "SKY")) {
     214        psError(PS_ERR_IO, false, "writing sky data\n");
     215        psFree(table);
     216        return false;
     217    }
     218
     219    psFree (table);
     220    psFree (header);
     221
     222    // second layer is the tangent plane
     223    header = psMetadataAlloc();
     224    psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "TANGENT_PLANE");
     225    psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "SKY");
     226    psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "RECTANGLE");
     227    psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "PROJECTION");
     228
     229    table = psArrayAllocEmpty (1);
     230    row = psMetadataAlloc ();
     231    psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "TANGENT_PLANE");
     232    psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "SKY");
     233
     234    psRegion *region = pmAstromFPAExtent (file->fpa);
     235    psMetadataAddF32(row,    PS_LIST_TAIL, "MINX",     PS_META_REPLACE, "range", region->x0);
     236    psMetadataAddF32(row,    PS_LIST_TAIL, "MAXX",     PS_META_REPLACE, "range", region->x1);
     237    psMetadataAddF32(row,    PS_LIST_TAIL, "MINY",     PS_META_REPLACE, "range", region->y0);
     238    psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
     239   
     240    psMetadataAddF32(row,    PS_LIST_TAIL, "XSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Xs * PM_DEG_RAD);
     241    psMetadataAddF32(row,    PS_LIST_TAIL, "YSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Ys * PM_DEG_RAD);
     242    psMetadataAddF32(row,    PS_LIST_TAIL, "XREF",     PS_META_REPLACE, "", file->fpa->toSky->R * PM_DEG_RAD);
     243    psMetadataAddF32(row,    PS_LIST_TAIL, "YREF",     PS_META_REPLACE, "", file->fpa->toSky->D * PM_DEG_RAD);
     244
     245    psArrayAdd (table, 100, row);
     246    psFree (row);
     247
     248    if (!psFitsWriteTable (file->fits, header, table, "TP")) {
     249        psError(PS_ERR_IO, false, "writing sky data\n");
     250        psFree(table);
     251        return false;
     252    }
     253
     254    psFree (table);
     255    psFree (header);
     256
     257    // third layer is the focal plane
     258    header = psMetadataAlloc();
     259    psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "FOCAL_PLANE");
     260    psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
     261    psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "RECTANGLE");
     262    psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "POLYNOMIAL");
     263
     264    table = psArrayAllocEmpty (1);
     265
     266    // XXX is this or the tpa region correct?
     267    region = pmAstromFPAExtent (file->fpa);
     268
     269    for (int i = 0; i <= file->fpa->toTPA->x->nX; i++) {
     270        for (int j = 0; j <= file->fpa->toTPA->x->nY; j++) {
     271            row = psMetadataAlloc ();
     272            psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "FOCAL_PLANE");
     273            psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
     274            psMetadataAddF32(row,    PS_LIST_TAIL, "MINX",     PS_META_REPLACE, "range", region->x0);
     275            psMetadataAddF32(row,    PS_LIST_TAIL, "MAXX",     PS_META_REPLACE, "range", region->x1);
     276            psMetadataAddF32(row,    PS_LIST_TAIL, "MINY",     PS_META_REPLACE, "range", region->y0);
     277            psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
     278   
     279            psMetadataAddS32(row,    PS_LIST_TAIL, "NX",       PS_META_REPLACE, "", file->fpa->toTPA->x->nX);
     280            psMetadataAddS32(row,    PS_LIST_TAIL, "NY",       PS_META_REPLACE, "", file->fpa->toTPA->x->nY);
     281            psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_X",   PS_META_REPLACE, "", file->fpa->toTPA->x->coeff[i][j]);
     282            psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_Y",   PS_META_REPLACE, "", file->fpa->toTPA->y->coeff[i][j]);
     283            psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_X",  PS_META_REPLACE, "", file->fpa->toTPA->x->coeffErr[i][j]);
     284            psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_Y",  PS_META_REPLACE, "", file->fpa->toTPA->y->coeffErr[i][j]);
     285
     286            psArrayAdd (table, 100, row);
     287            psFree (row);
     288        }
     289    }
     290
     291    if (!psFitsWriteTable (file->fits, header, table, "FP")) {
     292        psError(PS_ERR_IO, false, "writing sky data\n");
     293        psFree(table);
     294        return false;
     295    }
     296
     297    psFree (table);
     298    psFree (header);
     299
     300    // fourth layer holds the chips
     301    header = psMetadataAlloc();
     302    psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "FOCAL_PLANE");
     303    psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
     304    psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "RECTANGLE");
     305    psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "POLYNOMIAL");
     306
     307    table = psArrayAllocEmpty (1);
     308
     309    pmFPAview *view = pmFPAviewAlloc (0);
     310
     311    pmChip *chip = NULL;
     312    while ((chip = pmFPAviewNextChip (view, file->fpa, 1)) != NULL) {
     313
     314        psRegion *region = pmChipPixels (chip);
     315
     316        char *chiprule = psStringCopy ("CHIP.{CHIP.NAME}");
     317        char *chipname = pmFPAfileNameFromRule (chiprule, file, view);
     318
     319        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     320            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     321                row = psMetadataAlloc ();
     322                // xXXX set chip name based on metadata
     323                psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", chipname);
     324                psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
     325                psMetadataAddF32(row,    PS_LIST_TAIL, "MINX",     PS_META_REPLACE, "range", region->x0);
     326                psMetadataAddF32(row,    PS_LIST_TAIL, "MAXX",     PS_META_REPLACE, "range", region->x1);
     327                psMetadataAddF32(row,    PS_LIST_TAIL, "MINY",     PS_META_REPLACE, "range", region->y0);
     328                psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
     329   
     330                psMetadataAddS32(row,    PS_LIST_TAIL, "NX",       PS_META_REPLACE, "", chip->toFPA->x->nX);
     331                psMetadataAddS32(row,    PS_LIST_TAIL, "NY",       PS_META_REPLACE, "", chip->toFPA->x->nY);
     332                psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_X",   PS_META_REPLACE, "", chip->toFPA->x->coeff[i][j]);
     333                psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_Y",   PS_META_REPLACE, "", chip->toFPA->y->coeff[i][j]);
     334                psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_X",  PS_META_REPLACE, "", chip->toFPA->x->coeffErr[i][j]);
     335                psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_Y",  PS_META_REPLACE, "", chip->toFPA->y->coeffErr[i][j]);
     336                psArrayAdd (table, 100, row);
     337                psFree (row);
     338            }
     339        }
     340        psFree (chiprule);
     341        psFree (chipname);
     342    }
     343    psFree (view);
     344
     345    if (!psFitsWriteTable (file->fits, header, table, "CHIPS")) {
     346        psError(PS_ERR_IO, false, "writing sky data\n");
     347        psFree(table);
     348        return false;
     349    }
     350
     351    psFree (table);
     352    psFree (header);
    215353    return true;
    216354}
     
    267405
    268406    // EXTNAME for image header
    269     rule = psMetadataLookupStr(&status, menu, "ASTROMETRY");
     407    rule = psMetadataLookupStr(&status, menu, "REF.ASTROM");
    270408    if (!rule) {
    271         psError(PS_ERR_UNKNOWN, false, "missing entry for ASTROMETRY in EXTNAME.RULES in camera.config");
     409        psError(PS_ERR_UNKNOWN, false, "missing entry for REF.ASTROM in EXTNAME.RULES in camera.config");
    272410        return false;
    273411    }
  • trunk/psModules/src/camera/pmFPAfile.c

    r15629 r15665  
    439439    if (!strcasecmp (type, "HEADER"))     {
    440440        return PM_FPA_FILE_HEADER;
     441    }
     442    if (!strcasecmp (type, "ASTROM"))     {
     443        return PM_FPA_FILE_ASTROM;
    441444    }
    442445
     
    473476      case PM_FPA_FILE_HEADER:
    474477        return ("HEADER");
     478      case PM_FPA_FILE_ASTROM:
     479        return ("ASTROM");
    475480      default:
    476481        return ("NONE");
Note: See TracChangeset for help on using the changeset viewer.