IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7282


Ignore:
Timestamp:
Jun 1, 2006, 4:15:35 PM (20 years ago)
Author:
magnier
Message:

finish up the basic astrometry functions

Location:
trunk/psModules/src/astrom
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/pmAstrometryDistortion.c

    r7153 r7282  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-05-19 15:49:37 $
     9*  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-06-02 02:15:35 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    5757    // XXX we need to get the chip dimensions from the metadata somewhere
    5858    int Nx = 2;
    59     int Ny = 5;
     59    int Ny = 2;
    6060    int DX = 2080 / Nx;
    6161    int DY = 4600 / Ny;
    … …  
    108108                psVectorExtend (dP, 100, 1);
    109109                psVectorExtend (dQ, 100, 1);
    110             }
     110                Npts++;
     111            }
     112            if (Npts < 5)
     113                continue;
    111114
    112115            // fit the collection of positions and offsets with a local 1st order gradient
    … …  
    114117            psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    115118            psVector *mask = psVectorAlloc (Npts, PS_TYPE_MASK);
     119            mask->n = Npts;
     120            psVectorInit (mask, 0);
    116121
    117122            pmAstromGradient *grad = pmAstromGradientAlloc ();
    … …  
    141146bool pmAstromFitDistortion(pmFPA *fpa, psArray *grads, psMetadata *config)
    142147{
     148
     149    psPolynomial2D *localX = NULL;
     150    psPolynomial2D *localY = NULL;
    143151
    144152    // assign the gradient elements to psVectors for fitting
    … …  
    168176    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    169177    psVector *mask = psVectorAlloc (grads->n, PS_TYPE_MASK);
    170 
    171     // XXX the order of the gradient fits need to be 1 less than the distortion term
    172     // XXX determine the distortion order(s) from the fpa->toTP structure
    173 
    174     psPolynomial2D *localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
    175     psPolynomial2D *localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
    176     // XXX set masks based on fpa->toTP
    177 
     178    mask->n = grads->n;
     179    psVectorInit (mask, 0);
     180
     181    // the order of the gradient fits need to be 1 less than the distortion term
     182    // determine the gradient order(s) from the fpa->toTP structure
     183    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->x->nX-1, fpa->toTangentPlane->x->nY);
     184    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->x->nX,   fpa->toTangentPlane->x->nY-1);
     185
     186    // set masks based on fpa->toTP
     187    for (int i = 0; i < fpa->toTangentPlane->x->nX; i++) {
     188        for (int j = 0; j < fpa->toTangentPlane->x->nY; j++) {
     189            if (fpa->toTangentPlane->x->mask[i][j][0][0]) {
     190                localX->mask[i-1][j] = 1;
     191                localY->mask[i][j-1] = 1;
     192            }
     193        }
     194    }
     195
     196    // fit the local gradients in each direction
    178197    psVectorClipFitPolynomial2D (localX, stats, mask, 0xff, dPdL, NULL, L, M);
    179198    psVectorClipFitPolynomial2D (localY, stats, mask, 0xff, dPdM, NULL, L, M);
    180199
    181     // XXX update fpa->toTP distortion terms
    182 
     200    // XXX test: generate fit
     201    psVector *dPdLf = psPolynomial2DEvalVector (localX, L, M);
     202    psVector *dPdMf = psPolynomial2DEvalVector (localY, L, M);
     203
     204    // update fpa->toTP distortion terms
     205    fpa->toTangentPlane->x->coeff[0][0][0][0] = 0;
     206    for (int i = 1; i < fpa->toTangentPlane->x->nX; i++) {
     207        if (!fpa->toTangentPlane->x->mask[i][0][0][0]) {
     208            fpa->toTangentPlane->x->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
     209        }
     210    }
     211    for (int j = 1; j < fpa->toTangentPlane->x->nY; j++) {
     212        if (!fpa->toTangentPlane->x->mask[0][j][0][0]) {
     213            fpa->toTangentPlane->x->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
     214        }
     215    }
     216    for (int i = 1; i < fpa->toTangentPlane->x->nX; i++) {
     217        for (int j = 1; j < fpa->toTangentPlane->x->nY; j++) {
     218            if (!fpa->toTangentPlane->x->mask[i][j][0][0]) {
     219                fpa->toTangentPlane->x->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
     220            }
     221        }
     222    }
     223    psFree (localX);
     224    psFree (localY);
     225
     226    // the order of the gradient fits need to be 1 less than the distortion term
     227    // determine the gradient order(s) from the fpa->toTP structure
     228    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->y->nX-1, fpa->toTangentPlane->y->nY);
     229    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->y->nX,   fpa->toTangentPlane->y->nY-1);
     230
     231    // set masks based on fpa->toTP
     232    for (int i = 0; i < fpa->toTangentPlane->y->nX; i++) {
     233        for (int j = 0; j < fpa->toTangentPlane->y->nY; j++) {
     234            if (fpa->toTangentPlane->y->mask[i][j][0][0]) {
     235                localX->mask[i-1][j] = 1;
     236                localY->mask[i][j-1] = 1;
     237            }
     238        }
     239    }
     240
     241    // fit the local gradients in each direction
    183242    psVectorClipFitPolynomial2D (localX, stats, mask, 0xff, dQdL, NULL, L, M);
    184243    psVectorClipFitPolynomial2D (localY, stats, mask, 0xff, dQdM, NULL, L, M);
    185244
     245    // XXX test: generate fit
     246    psVector *dQdLf = psPolynomial2DEvalVector (localX, L, M);
     247    psVector *dQdMf = psPolynomial2DEvalVector (localY, L, M);
     248
     249    FILE *f = fopen ("grads.dat", "w");
     250    for (int i = 0; i < dPdLf->n; i++) {
     251        fprintf (f, "%f %f  | %f %f  %f %f | %f %f  %f %f\n",
     252                 L->data.F32[i], M->data.F32[i],
     253                 dPdL->data.F32[i], dPdM->data.F32[i],
     254                 dQdL->data.F32[i], dQdM->data.F32[i],
     255                 dPdLf->data.F32[i], dPdMf->data.F32[i],
     256                 dQdLf->data.F32[i], dQdMf->data.F32[i]);
     257    }
     258    fclose (f);
     259
     260    // update fpa->toTP distortion terms
     261    fpa->toTangentPlane->y->coeff[0][0][0][0] = 0;
     262    for (int i = 1; i < fpa->toTangentPlane->y->nX; i++) {
     263        if (!fpa->toTangentPlane->y->mask[i][0][0][0]) {
     264            fpa->toTangentPlane->y->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
     265        }
     266    }
     267    for (int j = 1; j < fpa->toTangentPlane->y->nY; j++) {
     268        if (!fpa->toTangentPlane->y->mask[0][j][0][0]) {
     269            fpa->toTangentPlane->y->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
     270        }
     271    }
     272    for (int i = 1; i < fpa->toTangentPlane->y->nX; i++) {
     273        for (int j = 1; j < fpa->toTangentPlane->y->nY; j++) {
     274            if (!fpa->toTangentPlane->y->mask[i][j][0][0]) {
     275                fpa->toTangentPlane->y->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
     276            }
     277        }
     278    }
     279    psFree (localX);
     280    psFree (localY);
     281
    186282    // XXX free unneeded structures
     283    psFree (dPdL);
     284    psFree (dPdM);
     285    psFree (dQdL);
     286    psFree (dQdM);
     287    psFree (L);
     288    psFree (M);
     289    psFree (stats);
     290    psFree (mask);
    187291
    188292    return true;

  • trunk/psModules/src/astrom/pmAstrometryObjects.c

    r7153 r7282  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-05-19 15:49:37 $
     10*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-06-02 02:15:35 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3434 ******************************************************************************/
    3535
     36// XXX an example (might be faster)
     37void pmAstromObjSortByChipX2 (psArray *objects)
     38{
     39
     40    int l,j,ir,i;
     41    float v1, v2;
     42    pmAstromObj *tmp;
     43    pmAstromObj *obj;
     44
     45    int N = objects->n;
     46
     47    if (N < 2)
     48        return;
     49    l = N >> 1;
     50    ir = N - 1;
     51    for (;;) {
     52        if (l > 0) {
     53            tmp = objects->data[--l];
     54        } else {
     55            tmp = objects->data[ir];
     56            objects->data[ir] = objects->data[0];
     57            if (--ir == 0) {
     58                objects->data[0] = tmp;
     59                return;
     60            }
     61        }
     62        i = l;
     63        j = (l << 1) + 1;
     64        while (j <= ir) {
     65            obj = objects->data[j];
     66            v1 = obj->chip->x;
     67            obj = objects->data[j+1];
     68            v2 = obj->chip->x;
     69            if (j < ir && v1 < v2)
     70                ++j;
     71            obj = objects->data[j];
     72            v1 = obj->chip->x;
     73            v2 = tmp->chip->x;
     74            if (v2 < v1) {
     75                objects->data[i] = objects->data[j];
     76                j += (i=j) + 1;
     77            } else
     78                j = ir + 1;
     79        }
     80        objects->data[i] = tmp;
     81    }
     82}
     83
    3684int pmAstromObjSortByFPX(
    3785    const void **a,
    … …  
    59107}
    60108
    61 
     109int pmAstromObjSortByChipX(
     110    const void **a,
     111    const void **b)
     112{
     113    if (PM_ASTROMETRYOBJECTS_DEBUG) {
     114        PS_ASSERT_PTR_NON_NULL(a, 0);
     115        PS_ASSERT_PTR_NON_NULL(*a, 0);
     116        PS_ASSERT_PTR_NON_NULL(b, 0);
     117        PS_ASSERT_PTR_NON_NULL(*b, 0);
     118    }
     119
     120    pmAstromObj *A = *(pmAstromObj **)a;
     121    pmAstromObj *B = *(pmAstromObj **)b;
     122
     123    psF32 diff = A->chip->x - B->chip->x;
     124    if (diff > FLT_EPSILON) {
     125        return (+1);
     126    }
     127
     128    if (diff < FLT_EPSILON) {
     129        return (-1);
     130    }
     131    return (0);
     132}
    62133
    63134/******************************************************************************
    … …  
    116187    st2 = psArraySort (st2, pmAstromObjSortByFPX);
    117188
     189    // pmAstromObjSortByChipX2 (st1);
     190    // pmAstromObjSortByChipX2 (st2);
     191
    118192    psArray *matches = psArrayAlloc (100);
    119193    matches->n = 0;
    … …  
    124198
    125199        dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
    126         if (dX < -RADIUS) {
     200        if (dX <= -RADIUS) {
    127201            i++;
    128202            continue;
    129203        }
    130         if (dX > +RADIUS) {
     204        if (dX >= +RADIUS) {
    131205            j++;
    132206            continue;
    … …  
    134208
    135209        jStart = j;
    136         while ((dX > -RADIUS) && (j < st2->n)) {
     210        while ((fabs(dX) < RADIUS) && (j < st2->n)) {
    137211
    138212            dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
    … …  
    155229        i++;
    156230    }
    157     psLogMsg (__func__, 3, "radius match: %d pairs\n", matches->n);
     231    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
     232    return (matches);
     233}
     234
     235
     236
     237/******************************************************************************
     238pmAstromRadiusMatchChip(st1, st2, config) : match objects on the chip
     239 ******************************************************************************/
     240psArray *pmAstromRadiusMatchChip(
     241    psArray *st1,
     242    psArray *st2,
     243    psMetadata *config)
     244{
     245    PS_ASSERT_PTR_NON_NULL(st1, NULL);
     246    PS_ASSERT_PTR_NON_NULL(st2, NULL);
     247    PS_ASSERT_PTR_NON_NULL(config, NULL);
     248    bool status;
     249    int jStart;
     250    double dX, dY, dR;
     251
     252    double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS");
     253    double RADIUS_SQR = PS_SQR (RADIUS);
     254
     255    // sort both lists by Focal Plane X coord
     256    st1 = psArraySort (st1, pmAstromObjSortByChipX);
     257    st2 = psArraySort (st2, pmAstromObjSortByChipX);
     258
     259    psArray *matches = psArrayAlloc (100);
     260    matches->n = 0;
     261
     262    int i = 0;
     263    int j = 0;
     264    while ((i < st1->n) && (j < st2->n)) {
     265
     266        dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
     267        if (dX <= -RADIUS) {
     268            i++;
     269            continue;
     270        }
     271        if (dX >= +RADIUS) {
     272            j++;
     273            continue;
     274        }
     275
     276        jStart = j;
     277        while ((dX > -RADIUS) && (j < st2->n)) {
     278
     279            dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
     280            dY = ((pmAstromObj *)st1->data[i])->chip->y - ((pmAstromObj *)st2->data[j])->chip->y;
     281            dR = dX*dX + dY*dY;
     282
     283            if (dR > RADIUS_SQR) {
     284                j++;
     285                continue;
     286            }
     287
     288            // got a match; add to output list
     289            pmAstromMatch *match = pmAstromMatchAlloc (i, j);
     290            psArrayAdd (matches, 100, match);
     291            psFree (match);
     292
     293            j++;
     294        }
     295        j = jStart;
     296        i++;
     297    }
     298    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
    158299    return (matches);
    159300}
    … …  
    235376    // fit chip-to-FPA transformation
    236377    psVectorClipFitPolynomial2D (map->x, stats, xMask, 0, x, wt, X, Y);
     378    psLogMsg ("psModules.astrom", 3, "x resid: %f +/- %f (%d of %d)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, x->n);
     379
    237380    psVectorClipFitPolynomial2D (map->y, stats, yMask, 0, y, wt, X, Y);
     381    psLogMsg ("psModules.astrom", 3, "y resid: %f +/- %f (%d of %d)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, y->n);
    238382
    239383    psFree (x);

  • trunk/psModules/src/astrom/pmAstrometryObjects.h

    r7153 r7282  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-05-19 15:49:37 $
     10*  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-06-02 02:15:35 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    9292 */
    9393psArray *pmAstromRadiusMatch(
     94    psArray *st1,
     95    psArray *st2,
     96    psMetadata *config
     97);
     98
     99psArray *pmAstromRadiusMatchChip(
    94100    psArray *st1,
    95101    psArray *st2,
Note: See TracChangeset for help on using the changeset viewer.