IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 30, 2006, 12:03:58 PM (20 years ago)
Author:
magnier
Message:

added pmAstromGridTweak, changed pmAstromObj elements to pointers

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r6872 r7005  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-04-17 18:01:04 $
     10*  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-04-30 22:03:58 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2020#include <math.h>
    2121#include "pslib.h"
     22#include "psVectorSmooth.h"
    2223#include "pmFPA.h"
    2324#include "pmAstrometryObjects.h"
    … …  
    4647    pmAstromObj *B = *(pmAstromObj **)b;
    4748
    48     psF32 diff = A->FP.x - B->FP.x;
     49    psF32 diff = A->FP->x - B->FP->x;
    4950    if (diff > FLT_EPSILON) {
    5051        return (+1);
    … …  
    121122    while ((i < st1->n) && (j < st2->n)) {
    122123
    123         dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x;
     124        dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
    124125        if (dX < -RADIUS) {
    125126            i++;
    … …  
    134135        while ((dX > -RADIUS) && (j < st2->n)) {
    135136
    136             dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x;
    137             dY = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y;
     137            dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
     138            dY = ((pmAstromObj *)st1->data[i])->FP->y - ((pmAstromObj *)st2->data[j])->FP->y;
    138139            dR = dX*dX + dY*dY;
    139140
    … …  
    146147            pmAstromMatch *match = pmAstromMatchAlloc (i, j);
    147148            psArrayAdd (matches, 100, match);
     149            psFree (match);
    148150
    149151            j++;
    … …  
    152154        i++;
    153155    }
     156    psLogMsg (__func__, 3, "radius match: %d pairs\n", matches->n);
    154157    return (matches);
    155158}
    … …  
    213216        refStar = ref->data[pair->i2];
    214217
    215         X->data.F32[i] = rawStar->chip.x;
    216         Y->data.F32[i] = rawStar->chip.y;
    217 
    218         x->data.F32[i] = refStar->FP.x;
    219         y->data.F32[i] = refStar->FP.y;
     218        X->data.F32[i] = rawStar->chip->x;
     219        Y->data.F32[i] = rawStar->chip->y;
     220
     221        x->data.F32[i] = refStar->FP->x;
     222        y->data.F32[i] = refStar->FP->y;
    220223
    221224        wt->data.F32[i] = 1.0;
    222225    }
    223226
    224     // no masking, no errors
    225     psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, wt, x, y);
    226     psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, wt, x, y);
     227    // constant errors
     228    psVector *xMask = psVectorAlloc (match->n, PS_TYPE_U8);
     229    psVector *yMask = psVectorAlloc (match->n, PS_TYPE_U8);
     230    xMask->n = yMask->n = match->n;
     231    psVectorInit (xMask, 0);
     232    psVectorInit (yMask, 0);
     233
     234    // fit chip-to-FPA transformation
     235    psVectorClipFitPolynomial2D (map->x, stats, xMask, 0, x, wt, X, Y);
     236    psVectorClipFitPolynomial2D (map->y, stats, yMask, 0, y, wt, X, Y);
     237
    227238    psFree (x);
    228239    psFree (y);
    … …  
    231242    psFree (wt);
    232243    psFree (stats);
     244    psFree (xMask);
     245    psFree (yMask);
    233246
    234247    return (map);
    235248}
    236 
    237249
    238250
    … …  
    264276        newObj = pmAstromObjCopy (oldObj);
    265277
    266         X = oldObj->FP.x - xCenter;
    267         Y = oldObj->FP.y - yCenter;
    268 
    269         newObj->FP.x = X*cs + Y*sn + xCenter;
    270         newObj->FP.y = Y*cs - X*sn + yCenter;
     278        X = oldObj->FP->x - xCenter;
     279        Y = oldObj->FP->y - yCenter;
     280
     281        newObj->FP->x = X*cs + Y*sn + xCenter;
     282        newObj->FP->y = Y*cs - X*sn + yCenter;
    271283
    272284        new->data[i] = newObj;
    … …  
    286298    }
    287299
     300    psFree(obj->pix);
     301    psFree(obj->cell);
     302    psFree(obj->chip);
     303    psFree(obj->FP);
     304    psFree(obj->TP);
     305    psFree(obj->sky);
     306
    288307    return;
    289308}
    … …  
    298317    psMemSetDeallocator (obj, (psFreeFunc) pmAstromObjFree);
    299318
    300     obj->pix.x = 0;
    301     obj->pix.y = 0;
    302 
    303     obj->FP.x = 0;
    304     obj->FP.y = 0;
    305 
    306     obj->TP.x = 0;
    307     obj->TP.y = 0;
    308 
    309     obj->sky.r = 0;
    310     obj->sky.d = 0;
    311 
    312     obj->Mag = 0;
     319    obj->pix  = psPlaneAlloc();
     320    obj->cell = psPlaneAlloc();
     321    obj->chip = psPlaneAlloc();
     322    obj->FP   = psPlaneAlloc();
     323    obj->TP   = psPlaneAlloc();
     324    obj->sky  = psSphereAlloc();
     325    obj->Mag  = 0;
    313326    obj->dMag = 0;
    314327
    … …  
    325338    pmAstromObj *obj = pmAstromObjAlloc();
    326339
    327     obj[0] = old[0];
     340    *obj->pix  = *old->pix;
     341    *obj->cell = *old->cell;
     342    *obj->chip = *old->chip;
     343    *obj->FP   = *old->FP;
     344    *obj->TP   = *old->TP;
     345    *obj->sky  = *old->sky;
    328346
    329347    return(obj);
    … …  
    450468        for (int j = 0; j < ref->n; j++) {
    451469            ob2 = (pmAstromObj *)ref->data[j];
    452             dX = ob1->FP.x - ob2->FP.x;
    453             dY = ob1->FP.y - ob2->FP.y;
    454 
    455             // 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);
     470            dX = ob1->FP->x - ob2->FP->x;
     471            dY = ob1->FP->y - ob2->FP->y;
     472
     473            // 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);
    456474            // find bin coordinates for this delta-delta
    457475            if (!AstromGridBin (&iX, &iY, dX, dY)) {
    … …  
    491509                    continue;
    492510
    493                 // this metric emphasize a narrow peak with lots of sources over one with few.
     511                // this metric emphasizes a narrow peak with lots of sources over one with few.
    494512                var = fabs((D2[j][i]/NP[j][i]) - PS_SQR(DX[j][i]/NP[j][i]) - PS_SQR(DY[j][i]/NP[j][i]));
    495                 metric = var / PS_SQR(PS_SQR(NP[j][i]));
    496 
    497                 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);
     513                metric = var / PS_SQR(NP[j][i]) / PS_SQR(NP[j][i]);
     514
     515                // 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);
    498516
    499517                if (metric < minMetric) {
    … …  
    513531        stats.nMatch    = NP[minY][minX];
    514532
     533        psFree (imStats);
    515534        // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better
    516535    }
     536    psFree (gridNP);
     537    psFree (gridDX);
     538    psFree (gridDY);
     539    psFree (gridD2);
    517540    return (stats);
    518541}
    … …  
    549572    for (int i = 0; i < raw->n; i++) {
    550573        obj = (pmAstromObj *)raw->data[i];
    551         xMin = PS_MIN (obj->FP.x, xMin);
    552         xMax = PS_MAX (obj->FP.x, xMax);
    553         yMin = PS_MIN (obj->FP.y, yMin);
    554         yMax = PS_MAX (obj->FP.y, yMax);
     574        xMin = PS_MIN (obj->FP->x, xMin);
     575        xMax = PS_MAX (obj->FP->x, xMax);
     576        yMin = PS_MIN (obj->FP->y, yMin);
     577        yMax = PS_MAX (obj->FP->y, yMax);
    555578    }
    556579    center.x = 0.5*(xMin + xMax);
    … …  
    565588        rot = pmAstromRotateObj (raw, center, angle);
    566589        newStat = pmAstromGridAngle (rot, ref, config);
     590        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));
     591
    567592        newStat.angle  = angle;
    568593        newStat.center = center;
    569594        if (newStat.minMetric < minStat.minMetric) {
    570595            minStat = newStat;
     596            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));
    571597        }
    572598        psFree (rot);
    573599    }
    574     fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric);
     600    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));
    575601    return (minStat);
    576602}
    577603
     604/******************************************************************************
     605pmAstromGridTweak(*raw, *ref, *recipe, stats): improve match for two star lists.
     606 ******************************************************************************/
     607pmAstromStats pmAstromGridTweak(
     608    psArray *raw,
     609    psArray *ref,
     610    psMetadata *recipe,
     611    pmAstromStats stats)
     612{
     613    bool status;
     614    pmAstromObj *ob1, *ob2;  // short-cut pointers to the objects
     615    double dX, dY;   // offset between a possible matched pair
     616    psArray *rot;
     617    int nBin, xBin, yBin;
     618
     619    rot = pmAstromRotateObj (raw, stats.center, stats.angle);
     620
     621    // sampling scale of the grid
     622    double tweakScale  = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.SCALE");
     623    double tweakRange  = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.RANGE");
     624    double tweakSmooth = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.SMOOTH");
     625    double tweakNsigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.NSIGMA");
     626
     627    nBin = 2*tweakRange / tweakScale;
     628    psVector *xHist = psVectorAlloc (nBin, PS_TYPE_F32);
     629    psVector *yHist = psVectorAlloc (nBin, PS_TYPE_F32);
     630    xHist->n = yHist->n = nBin;
     631    psVectorInit (xHist, 0);
     632    psVectorInit (yHist, 0);
     633
     634    // accumulate grids for focal plane (L,M) matches
     635    for (int i = 0; i < rot->n; i++) {
     636        ob1 = (pmAstromObj *)rot->data[i];
     637        for (int j = 0; j < ref->n; j++) {
     638            ob2 = (pmAstromObj *)ref->data[j];
     639            dX = ob1->FP->x - ob2->FP->x - stats.offset.x;
     640            dY = ob1->FP->y - ob2->FP->y - stats.offset.y;
     641
     642            xBin = (dX + tweakRange) / tweakScale;
     643            yBin = (dY + tweakRange) / tweakScale;
     644
     645            if (xBin < 0)
     646                continue;
     647            if (yBin < 0)
     648                continue;
     649            if (xBin >= nBin)
     650                continue;
     651            if (yBin >= nBin)
     652                continue;
     653
     654            xHist->data.F32[xBin] += 1.0;
     655            yHist->data.F32[yBin] += 1.0;
     656        }
     657    }
     658
     659    // smooth histgram vector with gaussian of 1sigma = radius
     660    psVectorSmooth (xHist, tweakSmooth, tweakNsigma);
     661    psVectorSmooth (yHist, tweakSmooth, tweakNsigma);
     662
     663    // select peak in x and in y
     664    xBin = yBin = 0;
     665    double xMax = 0;
     666    double yMax = 0;
     667    for (int i = 0; i < nBin; i++) {
     668        if (xHist->data.F32[i] > xMax) {
     669            xBin = i;
     670            xMax = xHist->data.F32[i];
     671        }
     672        if (yHist->data.F32[i] > yMax) {
     673            yBin = i;
     674            yMax = yHist->data.F32[i];
     675        }
     676    }
     677    double xPeak = xBin*tweakScale - tweakRange;
     678    double yPeak = yBin*tweakScale - tweakRange;
     679    psLogMsg (__func__, 3, "tweak peak by %f,%f\n", xPeak, yPeak);
     680
     681    // adjust offset by peak center
     682    pmAstromStats tweak = stats;
     683    tweak.offset.x += xPeak;
     684    tweak.offset.y += yPeak;
     685
     686    psFree (rot);
     687    psFree (xHist);
     688    psFree (yHist);
     689
     690    return tweak;
     691}
    578692
    579693/******************************************************************************
Note: See TracChangeset for help on using the changeset viewer.