IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:50:52 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/simtest_nebulous_branches
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/simtest_nebulous_branches

  • branches/simtest_nebulous_branches/psModules

  • branches/simtest_nebulous_branches/psModules/src/objects/pmSourceMatch.c

    r24262 r27840  
    88
    99#include "pmSource.h"
     10#include "pmErrorCodes.h"
    1011
    1112#include "pmSourceMatch.h"
    … …  
    112113    psFree(match->mag);
    113114    psFree(match->magErr);
     115    psFree(match->x);
     116    psFree(match->y);
    114117    psFree(match->image);
    115118    psFree(match->index);
     119    psFree(match->mask);
    116120}
    117121
    … …  
    124128    match->mag = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
    125129    match->magErr = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     130    match->x = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     131    match->y = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
    126132    match->image = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
    127133    match->index = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
    … …  
    133139void pmSourceMatchAdd(pmSourceMatch *match, // Match data
    134140                      float mag, float magErr, // Magnitude and error
     141                      float x, float y,        // Position
    135142                      int image, // Image index
    136143                      int index // Source index
    … …  
    141148    match->mag = psVectorExtend(match->mag, match->mag->nalloc, 1);
    142149    match->magErr = psVectorExtend(match->magErr, match->magErr->nalloc, 1);
     150    match->x = psVectorExtend(match->x, match->x->nalloc, 1);
     151    match->y = psVectorExtend(match->y, match->y->nalloc, 1);
    143152    match->image = psVectorExtend(match->image, match->image->nalloc, 1);
    144153    match->index = psVectorExtend(match->index, match->index->nalloc, 1);
    … …  
    147156    match->mag->data.F32[num] = mag;
    148157    match->magErr->data.F32[num] = magErr;
     158    match->x->data.F32[num] = x;
     159    match->y->data.F32[num] = y;
    149160    match->image->data.S32[num] = image;
    150161    match->index->data.S32[num] = index;
    … …  
    172183    for (int i = 0; i < numImages; i++) {
    173184        psArray *sources = sourceArrays->data[i]; // Sources in image
    174         if (!sources) {
     185        if (!sources || sources->n == 0) {
    175186            continue;
    176187        }
    … …  
    192203                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    193204                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    194                                  i, indices->data.S32[j]);
     205                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    195206                matches->data[j] = match;
    196207            }
    … …  
    212223                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    213224                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    214                                  i, indices->data.S32[j]);
     225                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    215226                matches->data[k] = match;
    216227            }
    … …  
    221232        } else {
    222233            // Match with the master list
    223             psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
     234            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, PS_TREE_EUCLIDEAN, xMaster, yMaster); // kd Tree
    224235            long numMatch = 0;          // Number of matches
    225236
    … …  
    238249                    pmSourceMatch *match = matches->data[index]; // Match data
    239250                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    240                                      i, indices->data.S32[j]);
     251                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    241252                    numMatch++;
    242253                } else {
    … …  
    244255                    pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    245256                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    246                                      i, indices->data.S32[j]);
     257                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    247258                    xMaster->data.F32[numMaster] = xImage->data.F32[j];
    248259                    yMaster->data.F32[numMaster] = yImage->data.F32[j];
    … …  
    262273        psFree(magImage);
    263274        psFree(magErrImage);
    264     }
     275        psFree(indices);
     276    }
     277
     278    psFree(xMaster);
     279    psFree(yMaster);
     280    psFree(boundsMaster);
    265281
    266282    if (cullSingles) {
    … …  
    304320        pmSourceMatch *match = matches->data[i]; // Match of interest
    305321        for (int j = 0; j < match->num && !source; j++) {
    306             if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) {
     322            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j]) ||
     323                !isfinite(match->x->data.F32[j]) || !isfinite(match->y->data.F32[j])) {
    307324                continue;
    308325            }
    … …  
    365382        double star = 0.0, starErr = 0.0; // Accumulators for star
    366383        for (int j = 0; j < match->num; j++) {
    367             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     384            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    368385                continue;
    369386            }
    … …  
    396413        pmSourceMatch *match = matches->data[i]; // Matched stars
    397414        for (int j = 0; j < match->num; j++) {
    398             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     415            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    399416                continue;
    400417            }
    … …  
    424441        pmSourceMatch *match = matches->data[i]; // Matched stars
    425442        for (int j = 0; j < match->num; j++) {
    426             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     443            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    427444                continue;
    428445            }
    … …  
    543560        pmSourceMatch *match = matches->data[i]; // Matched stars
    544561        for (int j = 0; j < match->num; j++) {
    545             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     562            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    546563                continue;
    547564            }
    … …  
    564581                if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
    565582                    numRejected++;
    566                     match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
     583                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_PHOT;
    567584                }
    568585            }
    … …  
    599616    int numImages = zp->n;              // Number of images
    600617    int numStars = matches->n;          // Number of stars
     618    psVector *badImage = psVectorAlloc(numImages, PS_TYPE_U8); // Bad image?
     619    psVectorInit(badImage, 0);
     620
     621    // Check for data integrity
     622    {
     623        psVector *num = psVectorAlloc(numImages, PS_TYPE_S32); // Number of stars per image
     624        psVectorInit(num, 0);
     625        for (int i = 0; i < numStars; i++) {
     626            pmSourceMatch *match = matches->data[i]; // Matched stars
     627            for (int j = 0; j < match->num; j++) {
     628                int index = match->image->data.U32[j]; // Image index
     629                psAssert(index >= 0 && index < numImages, "Bad index: %d", index);
     630                num->data.S32[index]++;
     631            }
     632        }
     633        int numGood = 0;                // Number of good images
     634        for (int i = 0; i < numImages; i++) {
     635            if (num->data.S32[i] == 0 || !isfinite(zp->data.F32[i])) {
     636                badImage->data.U8[i] = 0xFF;
     637                continue;
     638            }
     639            numGood++;
     640        }
     641        psFree(num);
     642        if (numGood == 0) {
     643            psError(PM_ERR_DATA, true, "No images with good stars.");
     644            psFree(badImage);
     645            return false;
     646        }
     647    }
     648
    601649    psVector *trans = psVectorAlloc(numImages, PS_TYPE_F32); // Transparencies for each image, magnitudes
    602650    psVectorInit(trans, 0.0);
    603651    psVector *photo = psVectorAlloc(numImages, PS_TYPE_U8); // Photometric determination for each image
    604652    psVectorInit(photo, 0);
    605     psVector *badImage = psVectorAlloc(numImages, PS_TYPE_U8); // Bad image?
    606     psVectorInit(badImage, 0);
    607653    psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
    608654
    … …  
    627673            return NULL;
    628674        }
    629         psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
     675        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric\n", numPhoto, numImages);
    630676
    631677        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1);
    632         psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
     678        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
    633679
    634680        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1);
    … …  
    649695            return NULL;
    650696        }
    651         psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages);
     697        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric\n", numPhoto, numImages);
    652698
    653699        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2);
    654         psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100);
     700        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
    655701
    656702        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2);
    … …  
    668714    return trans;
    669715}
     716
     717
     718// Iterate on the star positions and image shifts
     719// Returns the solution chi^2
     720static float sourceMatchRelastroIterate(psVector *xShift, psVector *yShift, // Shift for image
     721                                        psVector *xStar, psVector *yStar,   // Position for star
     722                                        const psArray *matches // Array of matches
     723    )
     724{
     725    psAssert(matches, "Need list of matches");
     726
     727    int numImages = xShift->n;          // Number of images
     728    int numStars = matches->n;          // Number of stars
     729
     730    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
     731             "Need shifts");
     732    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
     733    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
     734             "Need star positions");
     735    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
     736
     737    // Solve the star positions
     738    psVectorInit(xStar, NAN);
     739    psVectorInit(yStar, NAN);
     740    psVector *starMask = psVectorAlloc(numStars, PS_TYPE_U8); // Mask for stars
     741    psVectorInit(starMask, 0xFF);
     742    int numGoodStars = 0;               // Number of stars with good measurements
     743    for (int i = 0; i < numStars; i++) {
     744        pmSourceMatch *match = matches->data[i]; // Matched stars
     745        int numMeasurements = 0;        // Number of unmasked measurements for star
     746        double xSum = 0.0, ySum = 0.0;  // Accumulators for star
     747        for (int j = 0; j < match->num; j++) {
     748            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
     749                continue;
     750            }
     751            numMeasurements++;
     752            int index = match->image->data.U32[j]; // Image index
     753
     754            xSum += match->x->data.F32[j] - xShift->data.F32[index];
     755            ySum += match->y->data.F32[j] - yShift->data.F32[index];
     756        }
     757        if (numMeasurements > 1) {
     758            // It's only a good star (contributing to the chi^2) if there's more than 1 measurement
     759            numGoodStars++;
     760            xStar->data.F32[i] = xSum / numMeasurements;
     761            yStar->data.F32[i] = ySum / numMeasurements;
     762            starMask->data.U8[i] = 0;
     763        }
     764    }
     765
     766    // Solve for the shifts
     767    psVectorInit(xShift, 0.0);
     768    psVectorInit(yShift, 0.0);
     769    psVector *num = psVectorAlloc(numImages, PS_TYPE_S32);    // Number of stars
     770    psVectorInit(num, 0);
     771    for (int i = 0; i < numStars; i++) {
     772        if (starMask->data.U8[i]) {
     773            continue;
     774        }
     775        pmSourceMatch *match = matches->data[i]; // Matched stars
     776        for (int j = 0; j < match->num; j++) {
     777            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
     778                continue;
     779            }
     780            int index = match->image->data.U32[j]; // Image index
     781
     782            xShift->data.F32[index] += match->x->data.F32[j] - xStar->data.F32[i];
     783            yShift->data.F32[index] += match->y->data.F32[j] - yStar->data.F32[i];
     784            num->data.S32[index]++;
     785        }
     786    }
     787    for (int i = 0; i < numImages; i++) {
     788        xShift->data.F32[i] /= num->data.S32[i];
     789        yShift->data.F32[i] /= num->data.S32[i];
     790        psTrace("psModules.objects", 3, "Shift for image %d: %f,%f\n",
     791                i, xShift->data.F32[i], yShift->data.F32[i]);
     792    }
     793    psFree(num);
     794
     795    // Once more through to evaluate chi^2
     796    float chi2 = 0.0;                   // chi^2 for iteration
     797    int dof = 0;                        // Degrees of freedom
     798    for (int i = 0; i < numStars; i++) {
     799        pmSourceMatch *match = matches->data[i]; // Matched stars
     800        if (starMask->data.U8[i]) {
     801            continue;
     802        }
     803        for (int j = 0; j < match->num; j++) {
     804            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     805                continue;
     806            }
     807
     808            int index = match->image->data.U32[j]; // Image index
     809            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
     810            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
     811
     812            chi2 += PS_SQR(dx) + PS_SQR(dy);
     813            dof++;
     814        }
     815    }
     816    dof -= numGoodStars + numImages;
     817    chi2 /= dof;
     818
     819    return chi2;
     820}
     821
     822// Reject star measurements
     823// Returns the fraction of measurements that were rejected
     824static float sourceMatchRelastroReject(const psVector *xShift, const psVector *yShift, // Shifts for each image
     825                                       const psVector *xStar, const psVector *yStar, // Positions for each star
     826                                       const psArray *matches, // Array of matches
     827                                       float chi2,             // chi^2 from fit
     828                                       float rej               // Rejection threshold
     829                                )
     830{
     831    psAssert(matches, "Need list of matches");
     832
     833    int numImages = xShift->n;          // Number of images
     834    int numStars = matches->n;          // Number of stars
     835
     836    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
     837             "Need shifts");
     838    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
     839    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
     840             "Need star positions");
     841    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
     842
     843    int numRejected = 0;                // Number rejected
     844    int numMeasurements = 0;            // Number of measurements
     845
     846    float thresh = PS_SQR(rej) * chi2;    // Threshold for rejection
     847
     848    for (int i = 0; i < numStars; i++) {
     849        pmSourceMatch *match = matches->data[i]; // Matched stars
     850        for (int j = 0; j < match->num; j++) {
     851            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
     852                continue;
     853            }
     854            numMeasurements++;
     855            int index = match->image->data.U32[j]; // Image index
     856
     857            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
     858            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
     859
     860            if (PS_SQR(dx) + PS_SQR(dy) > thresh) {
     861                numRejected++;
     862                match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_ASTRO;
     863            }
     864        }
     865    }
     866
     867    return (float)numRejected / (float)numMeasurements;
     868}
     869
     870psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches
     871                               int numImages,          // Number of images
     872                               float tol, // Relative tolerance for convergence
     873                               int iter1, // Number of iterations for pass 1
     874                               float rej1, // Limit on rejection between iterations for pass 1
     875                               int iter2, // Number of iterations for pass 2
     876                               float rej2, // Limit on rejection between iterations for pass 2
     877                               float rejLimit // Limit on rejection between iterations
     878    )
     879{
     880    PS_ASSERT_ARRAY_NON_NULL(matches, NULL);
     881
     882    int numStars = matches->n;          // Number of stars
     883    psVector *xShift = psVectorAlloc(numImages, PS_TYPE_F32); // x shift for each image
     884    psVector *yShift = psVectorAlloc(numImages, PS_TYPE_F32); // y shift for each image
     885    psVectorInit(xShift, 0.0);
     886    psVectorInit(yShift, 0.0);
     887    psVector *xStar = psVectorAlloc(numStars, PS_TYPE_F32); // x position for each star
     888    psVector *yStar = psVectorAlloc(numStars, PS_TYPE_F32); // y position for each star
     889
     890    float chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); // chi^2 for solution
     891    psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
     892    float lastChi2 = INFINITY;          // chi^2 on last iteration
     893    float fracRej = INFINITY;           // Fraction of measurements rejected
     894
     895    // In the first passes, the shifts are not well deteremined: use high systematic error and
     896    // rejection thresholds
     897    for (int i = 0; i < iter1; i++) {
     898        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej1);
     899        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
     900
     901        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
     902        psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     903    }
     904
     905    for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
     906        lastChi2 = chi2;
     907
     908        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej2);
     909        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
     910
     911        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
     912        psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     913    }
     914
     915    psFree(xStar);
     916    psFree(yStar);
     917
     918    if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) {
     919        psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej);
     920    }
     921
     922    psArray *results = psArrayAlloc(numImages); // Array of results
     923    for (int i = 0; i < numImages; i++) {
     924        psVector *offset = results->data[i] = psVectorAlloc(2, PS_TYPE_F32); // Offset for image
     925        offset->data.F32[0] = xShift->data.F32[i];
     926        offset->data.F32[1] = yShift->data.F32[i];
     927    }
     928    psFree(xShift);
     929    psFree(yShift);
     930
     931    return results;
     932}
     933
Note: See TracChangeset for help on using the changeset viewer.