IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 18, 2006, 5:36:32 AM (20 years ago)
Author:
rhl
Message:

1/ Rewrote pmAstromRadiusMatch/pmAstromRadiusMatchChip so that they
do NOT modify their input lists (see bug 697). The two routines
now use a common worker routine. Tradeoffs: more memory allocation,
but no N ln N function calls inside qsort.
2/ Removed the functions that used to be passed to psArraySort
3/ Added const as appropriate, now that these routines don't
modify their input lists

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r7604 r7927  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-06-21 03:21:16 $
     10*  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-07-18 15:36:32 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2020#include <string.h>
    2121#include <math.h>
     22#include <assert.h>
    2223#include <unistd.h>   // for unlink
    2324#include <pslib.h>
    … …  
    2526#include "pmAstrometryObjects.h"
    2627
     28
    2729#define PM_ASTROMETRYOBJECTS_DEBUG 1
    2830/******************************************************************************
    29 pmAstromObjSortByFPX(**a, **b): sort by mag (descending)
     31pmAstromObjSortByMag(**a, **b): sort by mag (descending)
    3032 
    3133Is this a private routine?
    3234Should we do the early asserts?
    3335 ******************************************************************************/
    34 
    35 // XXX an example (might be faster)
    36 void pmAstromObjSortByChipX2 (psArray *objects)
    37 {
    38 
    39     int l,j,ir,i;
    40     float v1, v2;
    41     pmAstromObj *tmp;
    42     pmAstromObj *obj;
    43 
    44     int N = objects->n;
    45 
    46     if (N < 2)
    47         return;
    48     l = N >> 1;
    49     ir = N - 1;
    50     for (;;) {
    51         if (l > 0) {
    52             tmp = objects->data[--l];
    53         } else {
    54             tmp = objects->data[ir];
    55             objects->data[ir] = objects->data[0];
    56             if (--ir == 0) {
    57                 objects->data[0] = tmp;
    58                 return;
    59             }
    60         }
    61         i = l;
    62         j = (l << 1) + 1;
    63         while (j <= ir) {
    64             obj = objects->data[j];
    65             v1 = obj->chip->x;
    66             obj = objects->data[j+1];
    67             v2 = obj->chip->x;
    68             if (j < ir && v1 < v2)
    69                 ++j;
    70             obj = objects->data[j];
    71             v1 = obj->chip->x;
    72             v2 = tmp->chip->x;
    73             if (v2 < v1) {
    74                 objects->data[i] = objects->data[j];
    75                 j += (i=j) + 1;
    76             } else
    77                 j = ir + 1;
    78         }
    79         objects->data[i] = tmp;
    80     }
    81 }
    82 
    83 int pmAstromObjSortByFPX(
     36int pmAstromObjSortByMag(
    8437    const void **a,
    8538    const void **b)
    … …  
    9548    pmAstromObj *B = *(pmAstromObj **)b;
    9649
    97     psF32 diff = A->FP->x - B->FP->x;
    98     if (diff > FLT_EPSILON) {
    99         return (+1);
    100     }
    101 
    102     if (diff < FLT_EPSILON) {
    103         return (-1);
    104     }
    105     return (0);
    106 }
    107 
    108 int pmAstromObjSortByChipX(
    109     const void **a,
    110     const void **b)
    111 {
    112     if (PM_ASTROMETRYOBJECTS_DEBUG) {
    113         PS_ASSERT_PTR_NON_NULL(a, 0);
    114         PS_ASSERT_PTR_NON_NULL(*a, 0);
    115         PS_ASSERT_PTR_NON_NULL(b, 0);
    116         PS_ASSERT_PTR_NON_NULL(*b, 0);
    117     }
    118 
    119     pmAstromObj *A = *(pmAstromObj **)a;
    120     pmAstromObj *B = *(pmAstromObj **)b;
    121 
    122     psF32 diff = A->chip->x - B->chip->x;
    123     if (diff > FLT_EPSILON) {
    124         return (+1);
    125     }
    126 
    127     if (diff < FLT_EPSILON) {
    128         return (-1);
    129     }
    130     return (0);
    131 }
    132 
    133 /******************************************************************************
    134 pmAstromObjSortByMag(**a, **b): sort by mag (descending)
    135  
    136 Is this a private routine?
    137 Should we do the early asserts?
    138  ******************************************************************************/
    139 int pmAstromObjSortByMag(
    140     const void **a,
    141     const void **b)
    142 {
    143     if (PM_ASTROMETRYOBJECTS_DEBUG) {
    144         PS_ASSERT_PTR_NON_NULL(a, 0);
    145         PS_ASSERT_PTR_NON_NULL(*a, 0);
    146         PS_ASSERT_PTR_NON_NULL(b, 0);
    147         PS_ASSERT_PTR_NON_NULL(*b, 0);
    148     }
    149 
    150     pmAstromObj *A = *(pmAstromObj **)a;
    151     pmAstromObj *B = *(pmAstromObj **)b;
    152 
    15350    psF32 diff = A->Mag - B->Mag;
    15451    if (diff > FLT_EPSILON) {
    … …  
    16360}
    16461
    165 
    166 /******************************************************************************
    167 pmAstromRadiusMatch(st1, st2, config)
    168  ******************************************************************************/
    169 psArray *pmAstromRadiusMatch(
    170     psArray *st1,
    171     psArray *st2,
    172     psMetadata *config)
    173 {
    174     PS_ASSERT_PTR_NON_NULL(st1, NULL);
    175     PS_ASSERT_PTR_NON_NULL(st2, NULL);
    176     PS_ASSERT_PTR_NON_NULL(config, NULL);
    177     bool status;
     62/************************************************************************************************************/
     63/*
     64 * Working routine to match two lists (where x[12] are sorted), given psVectors of their coordinates and the
     65 * permutation used to sort in x
     66 */
     67static psArray *match_lists(const psVector *x1, const psVector *y1, // x/y coordinates of first set of objects
     68                            const psVector *x2, const psVector *y2, // x/y   "    "    "  second "   "  "   "
     69                            const psVector *sorted1, const psVector *sorted2, // mapping to original order
     70                            const double RADIUS) // matching radius
     71{
     72    psArray *matches = psArrayAlloc(x1->n);
     73    matches->n = 0;
     74
     75    const double RADIUS_SQR = PS_SQR(RADIUS);
     76    double dX, dY, dR;
     77
    17878    int jStart;
    179     double dX, dY, dR;
    180 
    181     double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS");
    182     double RADIUS_SQR = PS_SQR (RADIUS);
    183 
    184     // sort both lists by Focal Plane X coord
    185     st1 = psArraySort (st1, pmAstromObjSortByFPX);
    186     st2 = psArraySort (st2, pmAstromObjSortByFPX);
    187 
    188     // pmAstromObjSortByChipX2 (st1);
    189     // pmAstromObjSortByChipX2 (st2);
    190 
    191     psArray *matches = psArrayAlloc (100);
    192     matches->n = 0;
    193 
    194     int i = 0;
    195     int j = 0;
    196     while ((i < st1->n) && (j < st2->n)) {
    197 
    198         dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
     79    int i = 0, j = 0;
     80    while (i < x1->n && j < x2->n) {
     81        dX = x1->data.F64[i] - x2->data.F64[j];
    19982        if (dX <= -RADIUS) {
    20083            i++;
    … …  
    20790
    20891        jStart = j;
    209         while ((fabs(dX) < RADIUS) && (j < st2->n)) {
    210 
    211             dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
    212             dY = ((pmAstromObj *)st1->data[i])->FP->y - ((pmAstromObj *)st2->data[j])->FP->y;
     92        while (fabs(dX) < RADIUS && j < x2->n) {
     93
     94            dX = x1->data.F64[i] - x2->data.F64[j];
     95            dY = y1->data.F64[i] - y2->data.F64[j];
    21396            dR = dX*dX + dY*dY;
    21497
    … …  
    219102
    220103            // got a match; add to output list
    221             pmAstromMatch *match = pmAstromMatchAlloc (i, j);
     104            pmAstromMatch *match = pmAstromMatchAlloc (sorted1->data.S32[i], sorted2->data.S32[j]);
    222105            psArrayAdd (matches, 100, match);
    223106            psFree (match);
    … …  
    228111        i++;
    229112    }
    230     psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
     113
    231114    return (matches);
    232115}
    233 
    234 
    235 
    236 /******************************************************************************
    237 pmAstromRadiusMatchChip(st1, st2, config) : match objects on the chip
    238  ******************************************************************************/
    239 psArray *pmAstromRadiusMatchChip(
    240     psArray *st1,
    241     psArray *st2,
    242     psMetadata *config)
     116/******************************************************************************
     117pmAstromRadiusMatch(st1, st2, config)
     118 ******************************************************************************/
     119psArray *pmAstromRadiusMatch(
     120    const psArray *st1,
     121    const psArray *st2,
     122    const psMetadata *config)
    243123{
    244124    PS_ASSERT_PTR_NON_NULL(st1, NULL);
    245125    PS_ASSERT_PTR_NON_NULL(st2, NULL);
    246126    PS_ASSERT_PTR_NON_NULL(config, NULL);
    247     bool status;
    248     int jStart;
    249     double dX, dY, dR;
    250 
     127
     128    assert(st1->n == 0 || pmIsAstromObj(st1->data[0]));
     129    assert(st2->n == 0 || pmIsAstromObj(st2->data[0]));
     130
     131    bool status = false;
    251132    double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS");
    252     double RADIUS_SQR = PS_SQR (RADIUS);
    253 
    254     // sort both lists by Focal Plane X coord
    255     st1 = psArraySort (st1, pmAstromObjSortByChipX);
    256     st2 = psArraySort (st2, pmAstromObjSortByChipX);
    257 
    258     psArray *matches = psArrayAlloc (100);
    259     matches->n = 0;
    260 
    261     int i = 0;
    262     int j = 0;
    263     while ((i < st1->n) && (j < st2->n)) {
    264 
    265         dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
    266         if (dX <= -RADIUS) {
    267             i++;
    268             continue;
    269         }
    270         if (dX >= +RADIUS) {
    271             j++;
    272             continue;
    273         }
    274 
    275         jStart = j;
    276         while ((dX > -RADIUS) && (j < st2->n)) {
    277 
    278             dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
    279             dY = ((pmAstromObj *)st1->data[i])->chip->y - ((pmAstromObj *)st2->data[j])->chip->y;
    280             dR = dX*dX + dY*dY;
    281 
    282             if (dR > RADIUS_SQR) {
    283                 j++;
    284                 continue;
    285             }
    286 
    287             // got a match; add to output list
    288             pmAstromMatch *match = pmAstromMatchAlloc (i, j);
    289             psArrayAdd (matches, 100, match);
    290             psFree (match);
    291 
    292             j++;
    293         }
    294         j = jStart;
    295         i++;
    296     }
     133    if (!status) {
     134        psError(PS_ERR_IO, false, "Failed to lookup matching radius");
     135        return NULL;
     136    }
     137    /*
     138     * sort both lists by Focal Plane X coord; st1 first
     139     */
     140    psVector *x1 = psVectorAlloc(st1->n, PS_TYPE_F64);
     141    x1->n = x1->nalloc;
     142    for (int i = 0; i < st1->n; i++) {
     143        x1->data.F64[i] = ((pmAstromObj *)st1->data[i])->FP->x;
     144    }
     145    const psVector *sorted1 = psVectorSortIndex(NULL, x1);
     146    assert (sorted1->type.type == PS_TYPE_S32);
     147
     148    psVector *y1 = psVectorAlloc(st1->n, PS_TYPE_F64);
     149    y1->n = y1->nalloc;
     150    for (int i = 0; i < st1->n; i++) {
     151        x1->data.F64[i] = ((pmAstromObj *)st1->data[sorted1->data.S32[i]])->FP->x;
     152        y1->data.F64[i] = ((pmAstromObj *)st1->data[sorted1->data.S32[i]])->FP->y;
     153    }
     154
     155    // now st2
     156    psVector *x2 = psVectorAlloc(st2->n, PS_TYPE_F64);
     157    x2->n = x2->nalloc;
     158    for (int i = 0; i < st2->n; i++) {
     159        x2->data.F64[i] = ((pmAstromObj *)st2->data[i])->FP->x;
     160    }
     161    const psVector *sorted2 = psVectorSortIndex(NULL, x2);
     162
     163    psVector *y2 = psVectorAlloc(st2->n, PS_TYPE_F64);
     164    y2->n = y2->nalloc;
     165    for (int i = 0; i < st2->n; i++) {
     166        x2->data.F64[i] = ((pmAstromObj *)st2->data[sorted2->data.S32[i]])->FP->x;
     167        y2->data.F64[i] = ((pmAstromObj *)st2->data[sorted2->data.S32[i]])->FP->y;
     168    }
     169    /*
     170     * Do the work
     171     */
     172    psArray *matches = match_lists(x1, y1, x2, y2, sorted1, sorted2, RADIUS);
     173
     174    psFree(sorted1);
     175    psFree(sorted2);
     176    psFree(x1);
     177    psFree(y1);
     178    psFree(x2);
     179    psFree(y2);
     180
    297181    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
    298182    return (matches);
    299183}
    300184
    301 
     185/******************************************************************************
     186pmAstromRadiusMatchChip(st1, st2, config) : match objects on the chip
     187 ******************************************************************************/
     188psArray *pmAstromRadiusMatchChip(
     189    const psArray *st1,
     190    const psArray *st2,
     191    const psMetadata *config)
     192{
     193    PS_ASSERT_PTR_NON_NULL(st1, NULL);
     194    PS_ASSERT_PTR_NON_NULL(st2, NULL);
     195    PS_ASSERT_PTR_NON_NULL(config, NULL);
     196
     197    assert(st1->n == 0 || pmIsAstromObj(st1->data[0]));
     198    assert(st2->n == 0 || pmIsAstromObj(st2->data[0]));
     199
     200    bool status = false;
     201    double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS");
     202    if (!status) {
     203        psError(PS_ERR_IO, false, "Failed to lookup matching radius");
     204        return NULL;
     205    }
     206    /*
     207     * sort both lists by chip X coord; st1 first
     208     */
     209    psVector *x1 = psVectorAlloc(st1->n, PS_TYPE_F64);
     210    x1->n = x1->nalloc;
     211    for (int i = 0; i < st1->n; i++) {
     212        x1->data.F64[i] = ((pmAstromObj *)st1->data[i])->FP->x;
     213    }
     214    const psVector *sorted1 = psVectorSortIndex(NULL, x1);
     215    assert (sorted1->type.type == PS_TYPE_S32);
     216
     217    psVector *y1 = psVectorAlloc(st1->n, PS_TYPE_F64);
     218    y1->n = y1->nalloc;
     219    for (int i = 0; i < st1->n; i++) {
     220        x1->data.F64[i] = ((pmAstromObj *)st1->data[sorted1->data.S32[i]])->chip->x;
     221        y1->data.F64[i] = ((pmAstromObj *)st1->data[sorted1->data.S32[i]])->chip->y;
     222    }
     223
     224    // now st2
     225    psVector *x2 = psVectorAlloc(st2->n, PS_TYPE_F64);
     226    x2->n = x2->nalloc;
     227    for (int i = 0; i < st2->n; i++) {
     228        x2->data.F64[i] = ((pmAstromObj *)st2->data[i])->FP->x;
     229    }
     230    const psVector *sorted2 = psVectorSortIndex(NULL, x2);
     231
     232    psVector *y2 = psVectorAlloc(st2->n, PS_TYPE_F64);
     233    y2->n = y2->nalloc;
     234    for (int i = 0; i < st2->n; i++) {
     235        x2->data.F64[i] = ((pmAstromObj *)st2->data[sorted2->data.S32[i]])->chip->x;
     236        y2->data.F64[i] = ((pmAstromObj *)st2->data[sorted2->data.S32[i]])->chip->y;
     237    }
     238    /*
     239     * Do the work
     240     */
     241    psArray *matches = match_lists(x1, y1, x2, y2, sorted1, sorted2, RADIUS);
     242
     243    psFree(sorted1);
     244    psFree(sorted2);
     245    psFree(x1);
     246    psFree(y1);
     247    psFree(x2);
     248    psFree(y2);
     249
     250    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
     251    return (matches);
     252}
    302253
    303254/******************************************************************************
    … …  
    398349 ******************************************************************************/
    399350psArray *pmAstromRotateObj(
    400     psArray *old,
     351    const psArray *old,
    401352    psPlane center,
    402353    double angle)
    … …  
    406357    double X, Y;
    407358    pmAstromObj *newObj;
    408     pmAstromObj *oldObj;
     359    const pmAstromObj *oldObj;
    409360
    410361    psArray *new = psArrayAlloc (old->n);
    … …  
    434385
    435386/******************************************************************************
    436 pmAstromObjFree(obj)
    437  ******************************************************************************/
    438 static void pmAstromObjFree(pmAstromObj *obj)
     387astromObjFree(obj)
     388 ******************************************************************************/
     389static void astromObjFree(pmAstromObj *obj)
    439390{
    440391    if (obj == NULL) {
    … …  
    459410{
    460411    pmAstromObj *obj = psAlloc (sizeof(pmAstromObj));
    461     psMemSetDeallocator (obj, (psFreeFunc) pmAstromObjFree);
     412    psMemSetDeallocator (obj, (psFreeFunc)astromObjFree);
    462413
    463414    obj->pix  = psPlaneAlloc();
    … …  
    473424}
    474425
     426bool pmIsAstromObj(const psPtr ptr)
     427{
     428    return (psMemGetDeallocator(ptr) == (psFreeFunc)astromObjFree);
     429}
     430
     431
    475432
    476433/******************************************************************************
    477434pmAstromObjCopy(old)
    478435 ******************************************************************************/
    479 pmAstromObj *pmAstromObjCopy(pmAstromObj *old)
     436pmAstromObj *pmAstromObjCopy(const pmAstromObj *old)
    480437{
    481438    PS_ASSERT_PTR_NON_NULL(old, NULL);
    … …  
    557514 ******************************************************************************/
    558515pmAstromStats pmAstromGridAngle(
    559     psArray *raw,
    560     psArray *ref,
    561     psMetadata *config)
     516    const psArray *raw,
     517    const psArray *ref,
     518    const psMetadata *config)
    562519{
    563520    // XXX: What to do if input parameters are bad?
    … …  
    573530    int iX, iY;     // corresponding grid bin
    574531
    575     pmAstromObj *ob1, *ob2; // short-cut pointers to the objects
     532    const pmAstromObj *ob1, *ob2; // short-cut pointers to the objects
    576533    pmAstromStats stats;    // output match statistics
    577534
    … …  
    643600        // only check bins with at least 1/2 of max bin
    644601        int minNpts = 0.5*imStats->max;
    645         fprintf (stderr, "minNpts: %d, max: %d\n", minNpts, (int)(imStats->max));
     602        psTrace("psModule.astrom.grid.angle", 5, "minNpts: %d, max: %d", minNpts, (int)(imStats->max));
    646603
    647604        // find the 'best' bin
    … …  
    691648 ******************************************************************************/
    692649
    693 pmAstromStats pmAstromGridMatch(
    694     psArray *raw,
    695     psArray *ref,
    696     psMetadata *config)
     650pmAstromStats pmAstromGridMatch(const psArray *raw,
     651                                const psArray *ref,
     652                                const psMetadata *config)
    697653{
    698654    // XXX: What to do if input parameters are bad?
    … …  
    704660    bool status;
    705661    double xMin, xMax, yMin, yMax;
    706     pmAstromObj *obj;
     662    const pmAstromObj *obj;
    707663    psArray *rot;
    708664
    … …  
    732688        rot = pmAstromRotateObj (raw, center, angle);
    733689        newStat = pmAstromGridAngle (rot, ref, config);
    734         fprintf (stderr, "try %f : %f %f (%d pts, %f var, %f met (%f log))\n", angle, newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric, log10(newStat.minMetric));
     690        psTrace("psModule.astrom.grid.match", 5, "try %f : %f %f (%d pts, %f var, %f met (%f log))",
     691                angle, newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric,
     692                log10(newStat.minMetric));
    735693
    736694        newStat.angle  = angle;
    … …  
    738696        if (newStat.minMetric < minStat.minMetric) {
    739697            minStat = newStat;
    740             fprintf (stderr, "use %f : %f %f (%d pts, %f var, %f met (%f log))\n", angle, minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, minStat.minMetric, log10(minStat.minMetric));
     698            psTrace("psModule.astrom.grid.match", 5, "use %f : %f %f (%d pts, %f var, %f met (%f log))",
     699                    angle, minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, minStat.
     700                    minMetric, log10(minStat.minMetric));
    741701        }
    742702        psFree (rot);
    743703    }
    744     fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, log10(minStat.minMetric));
     704    psTrace("psModule.astrom.grid.match", 4, "best: %f %f (%d pts, %f var, %f met)",
     705            minStat.offset.x, minStat.offset.y, minStat.nMatch, minStat.minVar, log10(minStat.minMetric));
    745706    return (minStat);
    746707}
Note: See TracChangeset for help on using the changeset viewer.