IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 9, 2009, 2:53:12 PM (17 years ago)
Author:
Paul Price
Message:

Merging branches/pap (unconvolved stacks, reworked combinePixels function into clearer roles, only throw out most variant pixel on each iteration, throw suspect pixels out as first rejection step) stack development. I think I've got everything, but not entirely sure, since I've already merged this branch once before (for dual convolution).

Location:
trunk/psModules/src/objects
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects

  • trunk/psModules/src/objects/pmSourceMatch.c

    r25256 r26076  
    1717#define SOURCES_MAX_LEAF 2              // Maximum number of points on a tree leaf
    1818#define ARRAY_BUFFER 16                 // Buffer for array
     19
    1920
    2021//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    124125    match->mag = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
    125126    match->magErr = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     127    match->x = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     128    match->y = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
    126129    match->image = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
    127130    match->index = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
    … …  
    133136void pmSourceMatchAdd(pmSourceMatch *match, // Match data
    134137                      float mag, float magErr, // Magnitude and error
     138                      float x, float y,        // Position
    135139                      int image, // Image index
    136140                      int index // Source index
    … …  
    141145    match->mag = psVectorExtend(match->mag, match->mag->nalloc, 1);
    142146    match->magErr = psVectorExtend(match->magErr, match->magErr->nalloc, 1);
     147    match->x = psVectorExtend(match->x, match->x->nalloc, 1);
     148    match->y = psVectorExtend(match->y, match->y->nalloc, 1);
    143149    match->image = psVectorExtend(match->image, match->image->nalloc, 1);
    144150    match->index = psVectorExtend(match->index, match->index->nalloc, 1);
    … …  
    147153    match->mag->data.F32[num] = mag;
    148154    match->magErr->data.F32[num] = magErr;
     155    match->x->data.F32[num] = x;
     156    match->y->data.F32[num] = y;
    149157    match->image->data.S32[num] = image;
    150158    match->index->data.S32[num] = index;
    … …  
    192200                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    193201                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    194                                  i, indices->data.S32[j]);
     202                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    195203                matches->data[j] = match;
    196204            }
    … …  
    212220                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    213221                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    214                                  i, indices->data.S32[j]);
     222                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    215223                matches->data[k] = match;
    216224            }
    … …  
    238246                    pmSourceMatch *match = matches->data[index]; // Match data
    239247                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    240                                      i, indices->data.S32[j]);
     248                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    241249                    numMatch++;
    242250                } else {
    … …  
    244252                    pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    245253                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
    246                                      i, indices->data.S32[j]);
     254                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
    247255                    xMaster->data.F32[numMaster] = xImage->data.F32[j];
    248256                    yMaster->data.F32[numMaster] = yImage->data.F32[j];
    … …  
    304312        pmSourceMatch *match = matches->data[i]; // Match of interest
    305313        for (int j = 0; j < match->num && !source; j++) {
    306             if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) {
     314            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j]) ||
     315                !isfinite(match->x->data.F32[j]) || !isfinite(match->y->data.F32[j])) {
    307316                continue;
    308317            }
    … …  
    365374        double star = 0.0, starErr = 0.0; // Accumulators for star
    366375        for (int j = 0; j < match->num; j++) {
    367             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     376            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    368377                continue;
    369378            }
    … …  
    396405        pmSourceMatch *match = matches->data[i]; // Matched stars
    397406        for (int j = 0; j < match->num; j++) {
    398             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     407            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    399408                continue;
    400409            }
    … …  
    424433        pmSourceMatch *match = matches->data[i]; // Matched stars
    425434        for (int j = 0; j < match->num; j++) {
    426             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     435            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    427436                continue;
    428437            }
    … …  
    543552        pmSourceMatch *match = matches->data[i]; // Matched stars
    544553        for (int j = 0; j < match->num; j++) {
    545             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     554            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
    546555                continue;
    547556            }
    … …  
    564573                if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
    565574                    numRejected++;
    566                     match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
     575                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_PHOT;
    567576                }
    568577            }
    … …  
    627636            return NULL;
    628637        }
    629         psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
     638        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric\n", numPhoto, numImages);
    630639
    631640        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1);
    632         psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
     641        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
    633642
    634643        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1);
    … …  
    649658            return NULL;
    650659        }
    651         psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages);
     660        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric\n", numPhoto, numImages);
    652661
    653662        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2);
    654         psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100);
     663        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
    655664
    656665        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2);
    … …  
    668677    return trans;
    669678}
     679
     680
     681// Iterate on the star positions and image shifts
     682// Returns the solution chi^2
     683static float sourceMatchRelastroIterate(psVector *xShift, psVector *yShift, // Shift for image
     684                                        psVector *xStar, psVector *yStar,   // Position for star
     685                                        const psArray *matches // Array of matches
     686    )
     687{
     688    psAssert(matches, "Need list of matches");
     689
     690    int numImages = xShift->n;          // Number of images
     691    int numStars = matches->n;          // Number of stars
     692
     693    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
     694             "Need shifts");
     695    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
     696    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
     697             "Need star positions");
     698    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
     699
     700    // Solve the star positions
     701    psVectorInit(xStar, NAN);
     702    psVectorInit(yStar, NAN);
     703    psVector *starMask = psVectorAlloc(numStars, PS_TYPE_U8); // Mask for stars
     704    psVectorInit(starMask, 0xFF);
     705    int numGoodStars = 0;               // Number of stars with good measurements
     706    for (int i = 0; i < numStars; i++) {
     707        pmSourceMatch *match = matches->data[i]; // Matched stars
     708        int numMeasurements = 0;        // Number of unmasked measurements for star
     709        double xSum = 0.0, ySum = 0.0;  // Accumulators for star
     710        for (int j = 0; j < match->num; j++) {
     711            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
     712                continue;
     713            }
     714            numMeasurements++;
     715            int index = match->image->data.U32[j]; // Image index
     716
     717            xSum += match->x->data.F32[j] - xShift->data.F32[index];
     718            ySum += match->y->data.F32[j] - yShift->data.F32[index];
     719        }
     720        if (numMeasurements > 1) {
     721            // It's only a good star (contributing to the chi^2) if there's more than 1 measurement
     722            numGoodStars++;
     723            xStar->data.F32[i] = xSum / numMeasurements;
     724            yStar->data.F32[i] = ySum / numMeasurements;
     725            starMask->data.U8[i] = 0;
     726        }
     727    }
     728
     729    // Solve for the shifts
     730    psVectorInit(xShift, 0.0);
     731    psVectorInit(yShift, 0.0);
     732    psVector *num = psVectorAlloc(numImages, PS_TYPE_S32);    // Number of stars
     733    psVectorInit(num, 0);
     734    for (int i = 0; i < numStars; i++) {
     735        if (starMask->data.U8[i]) {
     736            continue;
     737        }
     738        pmSourceMatch *match = matches->data[i]; // Matched stars
     739        for (int j = 0; j < match->num; j++) {
     740            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
     741                continue;
     742            }
     743            int index = match->image->data.U32[j]; // Image index
     744
     745            xShift->data.F32[index] += match->x->data.F32[j] - xStar->data.F32[i];
     746            yShift->data.F32[index] += match->y->data.F32[j] - yStar->data.F32[i];
     747            num->data.S32[index]++;
     748        }
     749    }
     750    for (int i = 0; i < numImages; i++) {
     751        xShift->data.F32[i] /= num->data.S32[i];
     752        yShift->data.F32[i] /= num->data.S32[i];
     753        psTrace("psModules.objects", 3, "Shift for image %d: %f,%f\n",
     754                i, xShift->data.F32[i], yShift->data.F32[i]);
     755    }
     756    psFree(num);
     757
     758    // Once more through to evaluate chi^2
     759    float chi2 = 0.0;                   // chi^2 for iteration
     760    int dof = 0;                        // Degrees of freedom
     761    for (int i = 0; i < numStars; i++) {
     762        pmSourceMatch *match = matches->data[i]; // Matched stars
     763        if (starMask->data.U8[i]) {
     764            continue;
     765        }
     766        for (int j = 0; j < match->num; j++) {
     767            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     768                continue;
     769            }
     770
     771            int index = match->image->data.U32[j]; // Image index
     772            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
     773            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
     774
     775            chi2 += PS_SQR(dx) + PS_SQR(dy);
     776            dof++;
     777        }
     778    }
     779    dof -= numGoodStars + numImages;
     780    chi2 /= dof;
     781
     782    return chi2;
     783}
     784
     785// Reject star measurements
     786// Returns the fraction of measurements that were rejected
     787static float sourceMatchRelastroReject(const psVector *xShift, const psVector *yShift, // Shifts for each image
     788                                       const psVector *xStar, const psVector *yStar, // Positions for each star
     789                                       const psArray *matches, // Array of matches
     790                                       float chi2,             // chi^2 from fit
     791                                       float rej               // Rejection threshold
     792                                )
     793{
     794    psAssert(matches, "Need list of matches");
     795
     796    int numImages = xShift->n;          // Number of images
     797    int numStars = matches->n;          // Number of stars
     798
     799    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
     800             "Need shifts");
     801    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
     802    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
     803             "Need star positions");
     804    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
     805
     806    int numRejected = 0;                // Number rejected
     807    int numMeasurements = 0;            // Number of measurements
     808
     809    float thresh = PS_SQR(rej) * chi2;    // Threshold for rejection
     810
     811    for (int i = 0; i < numStars; i++) {
     812        pmSourceMatch *match = matches->data[i]; // Matched stars
     813        for (int j = 0; j < match->num; j++) {
     814            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
     815                continue;
     816            }
     817            numMeasurements++;
     818            int index = match->image->data.U32[j]; // Image index
     819
     820            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
     821            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
     822
     823            if (PS_SQR(dx) + PS_SQR(dy) > thresh) {
     824                numRejected++;
     825                match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_ASTRO;
     826            }
     827        }
     828    }
     829
     830    return (float)numRejected / (float)numMeasurements;
     831}
     832
     833psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches
     834                               int numImages,          // Number of images
     835                               float tol, // Relative tolerance for convergence
     836                               int iter1, // Number of iterations for pass 1
     837                               float rej1, // Limit on rejection between iterations for pass 1
     838                               int iter2, // Number of iterations for pass 2
     839                               float rej2, // Limit on rejection between iterations for pass 2
     840                               float rejLimit // Limit on rejection between iterations
     841    )
     842{
     843    PS_ASSERT_ARRAY_NON_NULL(matches, NULL);
     844
     845    int numStars = matches->n;          // Number of stars
     846    psVector *xShift = psVectorAlloc(numImages, PS_TYPE_F32); // x shift for each image
     847    psVector *yShift = psVectorAlloc(numImages, PS_TYPE_F32); // y shift for each image
     848    psVectorInit(xShift, 0.0);
     849    psVectorInit(yShift, 0.0);
     850    psVector *xStar = psVectorAlloc(numStars, PS_TYPE_F32); // x position for each star
     851    psVector *yStar = psVectorAlloc(numStars, PS_TYPE_F32); // y position for each star
     852
     853    float chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); // chi^2 for solution
     854    psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
     855    float lastChi2 = INFINITY;          // chi^2 on last iteration
     856    float fracRej = INFINITY;           // Fraction of measurements rejected
     857
     858    // In the first passes, the shifts are not well deteremined: use high systematic error and
     859    // rejection thresholds
     860    for (int i = 0; i < iter1; i++) {
     861        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej1);
     862        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
     863
     864        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
     865        psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     866    }
     867
     868    for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
     869        lastChi2 = chi2;
     870
     871        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej2);
     872        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
     873
     874        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
     875        psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     876    }
     877
     878    psFree(xStar);
     879    psFree(yStar);
     880
     881    if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) {
     882        psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej);
     883    }
     884
     885    psArray *results = psArrayAlloc(numImages); // Array of results
     886    for (int i = 0; i < numImages; i++) {
     887        psVector *offset = results->data[i] = psVectorAlloc(2, PS_TYPE_F32); // Offset for image
     888        offset->data.F32[0] = xShift->data.F32[i];
     889        offset->data.F32[1] = yShift->data.F32[i];
     890    }
     891    psFree(xShift);
     892    psFree(yShift);
     893
     894    return results;
     895}
     896
Note: See TracChangeset for help on using the changeset viewer.