IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 21, 2005, 11:03:53 AM (21 years ago)
Author:
eugene
Message:

various cleanups and fixes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/pmAstromGrid.c

    r5510 r5560  
    11# include "pmAstrom.h"
    22
    3 static double maxOffset;  // maximum allowed offset between lists, in raw pixels
     3static double maxOffpix;  // maximum allowed offset between lists, in raw pixels
    44static double Scale;      // grid pixel scale
    55static double Offset;     // deltas to pixels
    … …  
    77// local function to convert x,y coords to grid bins
    88// it requires the globals defined above
    9 static bool AstromGridBin (int *pOut, int *qOut, double dP, int dQ) {
    10 
    11     if (fabs(dP) > maxOffset) return false;
    12     if (fabs(dQ) > maxOffset) return false;
    13 
    14     *pOut = dP / Scale + Offset;
    15     *qOut = dQ / Scale + Offset;
     9static bool AstromGridBin (int *dx, int *dy, double dX, double dY) {
     10
     11    if (fabs(dX) > maxOffpix) return false;
     12    if (fabs(dY) > maxOffpix) return false;
     13
     14    *dx = dX / Scale + Offset;
     15    *dy = dY / Scale + Offset;
    1616    return true;
    1717}
    1818
    19 // match two star lists
    20 pmAstromGridMatchStats pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config) {
    21 
    22     double xMin, xMax, yMin, yMax;
    23     pmAstromObj *ob1, *ob2;
    24 
    25     pmAstromGridMatchStat minStat, newStat;
    26     psPlane center;
    27 
    28     // find center of the st2 field (focal-plane coords)
    29     xMin = yMin = +1e10;
    30     xMax = yMax = -1e10;
    31     for (int i = 0; i < st2->n; i++) {
    32         ob2 = (pmAstromObj *)st2->data[i];
    33         xMin = PS_MIN (ob2->FP.x, xMin);
    34         xMax = PS_MAX (ob2->FP.x, xMax);
    35         yMin = PS_MIN (ob2->FP.y, yMin);
    36         yMax = PS_MAX (ob2->FP.y, yMax);
    37     }
    38     center.x = 0.5*(xMin + xMax);
    39     center.y = 0.5*(yMin + yMax);
    40 
    41     double minAngle = psMetadataLookupF32 (&status, config, "GRID.MIN.ANGLE");
    42     double maxAngle = psMetadataLookupF32 (&status, config, "GRID.MAX.ANGLE");
    43     double delAngle = psMetadataLookupF32 (&status, config, "GRID.DEL.ANGLE");
    44 
    45     minStat.minMetric = 1e10;
    46     for (angle = minAngle; angle <= maxAngle; angle += delAngle) {
    47         st2r = pmAstromRotateObj (st2, center, angle);
    48         newStat = pmAstromGridMatchAngle (st1, st2r, config);
    49         newStat.angle  = angle;
    50         newStat.center = center;
    51         if (newStat.minMetric < minStat.minMetric) {
    52             minStat = newStat;
    53         }
    54         psFree (st2r);
    55     }
    56     return (minStat);
    57 }
    58 
    5919// match the two lists using the binned delta-delta max
    60 pmAstromGridMatchStats pmAstromGridMatchAngle (psArray *st1, psArray *st2, psMetadata *config) {
    61 
     20pmAstromStats pmAstromGridAngle (psArray *raw, psArray *ref, psMetadata *config) {
     21
     22    bool status;
    6223    int nPix;       // size of matching grid
     24    int nPixHalf;   // half-size of matching grid
    6325    double dX, dY;  // offset between a possible matched pair
    6426    int iX, iY;     // corresponding grid bin
    6527
    66     pmAstromObj *ob1, *ob2 // short-cut pointers to the objects
    67     pmAstromGridMatchStats matchStats; // output match statistics
     28    pmAstromObj *ob1, *ob2; // short-cut pointers to the objects
     29    pmAstromStats stats;    // output match statistics
    6830
    6931    // max allowed offset in either X or Y directions
    70     double gridOffset = psMetadataLookupF32 (&status, config, "GRID.OFFSET");
     32    double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET");
    7133
    7234    // sampling scale of the grid
    73     double gridScale  = psMetadataLookupF32 (&status, config, "GRID.SCALE");
     35    double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
    7436
    7537    // set the static scaling factors
    76     nPix = (int)(gridOffset / gridScale + 0.5);  // half-grid
    77     nPix = 2*nPix + 1;
     38    nPixHalf = (int)(gridOffset / gridScale + 0.5);  // half-grid
     39    nPix = 2*nPixHalf + 1;                           // full grid width
    7840
    7941    // these are globals used by p_pmAstromGridBin
    80     maxOffpix = gridScale * (nPix + 0.5);
     42    maxOffpix = gridScale * (nPixHalf + 0.5);            // max offset from true center
    8143    Offset    = maxOffpix / gridScale;
    8244    Scale     = gridScale;
    8345
    84     // XXX EAM : can we assume the allocated image is init-ed?
    85     psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_S32);
     46    // images used as accumulators for the loop below
     47    psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_U32);
    8648    psImage *gridDX = psImageAlloc (nPix, nPix, PS_TYPE_F32);
    8749    psImage *gridDY = psImageAlloc (nPix, nPix, PS_TYPE_F32);
    8850    psImage *gridD2 = psImageAlloc (nPix, nPix, PS_TYPE_F32);
    89     psImageInit (gridNP);
    90     psImageInit (gridDX);
    91     psImageInit (gridDY);
    92     psImageInit (gridD2);
     51    psImageInit (gridNP, 0);
     52    psImageInit (gridDX, 0);
     53    psImageInit (gridDY, 0);
     54    psImageInit (gridD2, 0);
    9355
    9456    // short-cut names for grid images
    95     psS32 *NP = gridNP->data.S32;
    96     psF32 *DX = gridDP->data.F32;
    97     psF32 *DY = gridDQ->data.F32;
    98     psF32 *D2 = gridD2->data.F32;
     57    psS32 **NP = gridNP->data.S32;
     58    psF32 **DX = gridDX->data.F32;
     59    psF32 **DY = gridDY->data.F32;
     60    psF32 **D2 = gridD2->data.F32;
    9961
    10062    // accumulate grids for focal plane (L,M) matches
    101     for (int i = 0; i < st1->n; i++) {
    102         ob1 = (pmAstromObj *)st1->data[i];
    103         for (int j = 0; j < st2->n; j++) {
    104             ob2 = (pmAstromObj *)st2->data[i];
     63    for (int i = 0; i < raw->n; i++) {
     64        ob1 = (pmAstromObj *)raw->data[i];
     65        for (int j = 0; j < ref->n; j++) {
     66            ob2 = (pmAstromObj *)ref->data[j];
    10567            dX = ob1->FP.x - ob2->FP.x;
    10668            dY = ob1->FP.y - ob2->FP.y;
    10769
     70            // fprintf (stderr, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP.x, ob2->FP.x, ob1->FP.y, ob2->FP.y);
    10871            // find bin coordinates for this delta-delta
    109             if (!p_pmAstromGridBin (&iX, &iY, dX, dY)) {
     72            if (!AstromGridBin (&iX, &iY, dX, dY)) {
    11073                continue; // matched pair is too far offset
    111 
    11274            }
    11375
    11476            // accumulate bin stats
    11577            NP[iY][iX] ++;
    116             DP[iY][iX] += dX;
    117             DQ[iY][iX] += dY;
     78            DX[iY][iX] += dX;
     79            DY[iY][iX] += dY;
    11880            D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY);
    11981        }
    … …  
    12486        double minMetric = 1e10;
    12587        double minVar = 1e10;
    126         double minX = -1;
    127         double minY = -1;
     88        int minX = -1;
     89        int minY = -1;
    12890        double metric, var;
    12991
    13092        // find the max pixel
    131         psStats *stats = psStatsAlloc (PS_STAT_MAX);
    132         stats = psImageStats (stats, gridNP);
     93        psStats *imStats = psStatsAlloc (PS_STAT_MAX);
     94        imStats = psImageStats (imStats, gridNP, NULL, 0);
    13395
    13496        // only check bins with at least 1/2 of max bin
    135         int minNpt = 0.5*stats->max;
     97        int minNpts = 0.5*imStats->max;
     98        fprintf (stderr, "minNpts: %d, max: %d\n", minNpts, (int)(imStats->max));
    13699
    137100        // find the 'best' bin
    138         for (int j = 0; j < gridNP->nY; j++) {
    139             for (int i = 0; i < gridNP->nX; i++) {
     101        for (int j = 0; j < gridNP->numRows; j++) {
     102            for (int i = 0; i < gridNP->numCols; i++) {
    140103
    141104                if (NP[j][i] < minNpts) continue;
    … …  
    145108                metric = var / PS_SQR(PS_SQR(NP[j][i]));
    146109           
     110                fprintf (stderr, "try : %f %f (%d pts, %f var, %f met)\n", DX[j][i]/NP[j][i], DY[j][i]/NP[j][i], NP[j][i], var, metric);
     111
    147112                if (metric < minMetric) {
    148113                    minMetric = metric;
    … …  
    155120
    156121        // convert the bin to delta-delta
    157         matchStats.offset.x  = DP[minY][minX] / NP[minY][minX];
    158         matchStats.offset.y  = DQ[minY][minX] / NP[minY][minX];
    159         matchStats.minMetric = minMetric;
    160         matchStats.minVar    = minVar;
    161         matchStats.nMatch    = NP[minY][minX];
    162     }
    163     return (matchStats);
     122        stats.offset.x  = DX[minY][minX] / NP[minY][minX];
     123        stats.offset.y  = DY[minY][minX] / NP[minY][minX];
     124        stats.minMetric = minMetric;
     125        stats.minVar    = minVar;
     126        stats.nMatch    = NP[minY][minX];
     127
     128        // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better
     129    }
     130    return (stats);
     131}
     132
     133// match two star lists
     134pmAstromStats pmAstromGridMatch (psArray *raw, psArray *ref, psMetadata *config) {
     135
     136    bool status;
     137    double xMin, xMax, yMin, yMax;
     138    pmAstromObj *obj;
     139    psArray *rot;
     140
     141    pmAstromStats minStat, newStat;
     142    psPlane center;
     143
     144    // find center of the raw field (focal-plane coords)
     145    xMin = yMin = +1e10;
     146    xMax = yMax = -1e10;
     147    for (int i = 0; i < raw->n; i++) {
     148        obj = (pmAstromObj *)raw->data[i];
     149        xMin = PS_MIN (obj->FP.x, xMin);
     150        xMax = PS_MAX (obj->FP.x, xMax);
     151        yMin = PS_MIN (obj->FP.y, yMin);
     152        yMax = PS_MAX (obj->FP.y, yMax);
     153    }
     154    center.x = 0.5*(xMin + xMax);
     155    center.y = 0.5*(yMin + yMax);
     156
     157    double minAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE");
     158    double maxAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE");
     159    double delAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE");
     160
     161    minStat.minMetric = 1e10;
     162    for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
     163        rot = pmAstromRotateObj (raw, center, angle);
     164        newStat = pmAstromGridAngle (rot, ref, config);
     165        newStat.angle  = angle;
     166        newStat.center = center;
     167        if (newStat.minMetric < minStat.minMetric) {
     168            minStat = newStat;
     169        }
     170        psFree (rot);
     171    }
     172    fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric);
     173    return (minStat);
    164174}
    165175
    166176// apply the measured FPA offset and rotation (stat) to the fpa astrom structures
    167 psFPA *pmAstromGridApply (psPlaneTransform *map, pmAstromGridMatchStat stat) {
    168 
    169     // stat.angle, stat.center, stat.offse
    170     // I think i need to know the center reference....
    171 
    172     return (fpa);
     177psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat) {
     178
     179    double cs = cos (stat.angle);
     180    double sn = sin (stat.angle);
     181   
     182    double dx = (map->x->coeff[0][0] - stat.center.x);
     183    double dy = (map->y->coeff[0][0] - stat.center.y);
     184
     185    // new offset
     186    map->x->coeff[0][0] =  cs*dx + sn*dy - stat.offset.x + stat.center.x;
     187    map->y->coeff[0][0] = -sn*dx + cs*dy - stat.offset.y + stat.center.y;
     188
     189    // original rotation matrix
     190    double pc1_1 = map->x->coeff[1][0];
     191    double pc1_2 = map->x->coeff[0][1];
     192    double pc2_1 = map->y->coeff[1][0];
     193    double pc2_2 = map->y->coeff[0][1];
     194
     195    // new rotation matrix
     196    map->x->coeff[1][0] = +cs*pc1_1 + sn*pc2_1;
     197    map->x->coeff[0][1] = +cs*pc1_2 + sn*pc2_2;
     198    map->y->coeff[1][0] = -sn*pc1_1 + cs*pc2_1;
     199    map->y->coeff[0][1] = -sn*pc1_2 + cs*pc2_2;
     200
     201    return (map);
    173202}
    174203
Note: See TracChangeset for help on using the changeset viewer.