IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 17, 2009, 12:07:42 PM (17 years ago)
Author:
beaumont
Message:

merged with head

Location:
branches/cnb_branches/cnb_branch_20090301/psModules
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/cnb_branches/cnb_branch_20090301/psModules

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmSourceMatch.c

    r21266 r23351  
    1616#define SOURCE_FAINTEST 50.0            // Faintest magnitude to consider
    1717#define SOURCES_MAX_LEAF 2              // Maximum number of points on a tree leaf
    18 
     18#define ARRAY_BUFFER 16                 // Buffer for array
    1919
    2020//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     
    6363        pmSource *source = sources->data[i]; // Source of interest
    6464        if (!source) continue;
    65         if (source->mode & SOURCE_MASK) continue;
    66         if (!isfinite(source->psfMag)) continue;
    67         if (!isfinite(source->errMag)) continue;
    68         if (source->psfMag > SOURCE_FAINTEST) continue;
     65        if (source->mode & SOURCE_MASK) continue;
     66        if (!isfinite(source->psfMag)) continue;
     67        if (!isfinite(source->errMag)) continue;
     68        if (source->psfMag > SOURCE_FAINTEST) continue;
    6969
    7070        float xSrc, ySrc;               // Coordinates of source
     
    113113}
    114114
    115 pmSourceMatch *pmSourceMatchAlloc(int num // Maximum number of images
    116     )
     115pmSourceMatch *pmSourceMatchAlloc(void)
    117116{
    118117    pmSourceMatch *match = psAlloc(sizeof(pmSourceMatch)); // Match data
     
    120119
    121120    match->num = 0;
    122     match->mag = psVectorAllocEmpty(num, PS_TYPE_F32);
    123     match->magErr = psVectorAllocEmpty(num, PS_TYPE_F32);
    124     match->image = psVectorAllocEmpty(num, PS_TYPE_U32);
    125     match->index = psVectorAllocEmpty(num, PS_TYPE_U32);
    126     match->mask = psVectorAllocEmpty(num, PS_TYPE_VECTOR_MASK);
     121    match->mag = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     122    match->magErr = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     123    match->image = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
     124    match->index = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
     125    match->mask = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_VECTOR_MASK);
    127126
    128127    return match;
     
    136135{
    137136    int num = match->num;               // Number of matches
    138     psAssert(num <= match->mag->nalloc, "Too many matches.");
    139     psAssert(num <= match->magErr->nalloc, "Too many matches.");
    140     psAssert(num <= match->image->nalloc, "Too many matches.");
    141     psAssert(num <= match->index->nalloc, "Too many matches.");
     137
     138    match->mag = psVectorExtend(match->mag, match->mag->nalloc, 1);
     139    match->magErr = psVectorExtend(match->magErr, match->magErr->nalloc, 1);
     140    match->image = psVectorExtend(match->image, match->image->nalloc, 1);
     141    match->index = psVectorExtend(match->index, match->index->nalloc, 1);
     142    match->mask = psVectorExtend(match->mask, match->mask->nalloc, 1);
    142143
    143144    match->mag->data.F32[num] = mag;
     
    148149    match->num++;
    149150
    150     match->mag->n = match->magErr->n = match->image->n = match->index->n = match->mask->n = match->num;
    151 }
    152 
    153 
    154 psArray *pmSourceMatchSources(const psArray *sourceArrays, float radius)
     151    return;
     152}
     153
     154
     155psArray *pmSourceMatchSources(const psArray *sourceArrays, float radius, bool cullSingles)
    155156{
    156157    PS_ASSERT_ARRAY_NON_NULL(sourceArrays, NULL);
     
    185186            matches = psArrayAlloc(numSources);
    186187            for (int j = 0; j < numSources; j++) {
    187                 pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data
     188                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    188189                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
    189190                matches->data[j] = match;
    190191            }
    191192            numMaster += numSources;
    192             goto match_image_done;
    193         }
    194 
    195         if (boundsImage->x0 > boundsMaster->x1 || boundsImage->x1 < boundsMaster->x0 ||
    196             boundsImage->y0 > boundsMaster->y1 || boundsImage->y1 < boundsMaster->y0) {
     193        } else if (boundsImage->x0 > boundsMaster->x1 || boundsImage->x1 < boundsMaster->x0 ||
     194                   boundsImage->y0 > boundsMaster->y1 || boundsImage->y1 < boundsMaster->y0) {
     195            // Bounds don't overlap --- can just add everything in to the master list
    197196            psTrace("psModules.objects", 7, "Bounds don't overlap\n");
    198             // Add everything in to the master list
    199197            long size = numMaster + numSources; // New size
    200198            xMaster = psVectorRealloc(xMaster, size);
     
    207205                   numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
    208206            for (int j = 0, k = numMaster; j < numSources; j++, k++) {
    209                 pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data
     207                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    210208                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
    211209                matches->data[k] = match;
     
    215213            yMaster->n = size;
    216214            matches->n = size;
    217 
    218             goto match_image_done;
    219         }
    220 
    221         // Match with the master list
    222 
    223         psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
    224         long numMatch = 0;              // Number of matches
    225 
    226         long size = numMaster + numSources; // New size
    227         xMaster = psVectorRealloc(xMaster, size);
    228         yMaster = psVectorRealloc(yMaster, size);
    229         matches = psArrayRealloc(matches, size);
    230 
    231         psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates for tree lookup
    232         for (int j = 0; j < numSources; j++) {
    233             coords->data.F32[0] = xImage->data.F32[j];
    234             coords->data.F32[1] = yImage->data.F32[j];
    235             long index = psTreeNearestWithin(tree, coords, radius); // Match index
    236             if (index >= 0) {
    237                 // Record the match
    238                 pmSourceMatch *match = matches->data[index]; // Match data
    239                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
    240                 numMatch++;
    241             } else {
    242                 // Add to the master list
    243                 pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data
    244                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
    245                 xMaster->data.F32[numMaster] = xImage->data.F32[j];
    246                 yMaster->data.F32[numMaster] = yImage->data.F32[j];
    247                 matches->data[numMaster] = match;
    248                 numMaster++;
    249                 xMaster->n = yMaster->n = matches->n = numMaster;
    250             }
    251 
    252         }
    253         psFree(coords);
    254         psFree(tree);
    255 
    256 
    257     match_image_done:
     215        } else {
     216            // Match with the master list
     217            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
     218            long numMatch = 0;          // Number of matches
     219
     220            long size = numMaster + numSources; // New size
     221            xMaster = psVectorRealloc(xMaster, size);
     222            yMaster = psVectorRealloc(yMaster, size);
     223            matches = psArrayRealloc(matches, size);
     224
     225            psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates for tree lookup
     226            for (int j = 0; j < numSources; j++) {
     227                coords->data.F32[0] = xImage->data.F32[j];
     228                coords->data.F32[1] = yImage->data.F32[j];
     229                long index = psTreeNearestWithin(tree, coords, radius); // Match index
     230                if (index >= 0) {
     231                    // Record the match
     232                    pmSourceMatch *match = matches->data[index]; // Match data
     233                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     234                    numMatch++;
     235                } else {
     236                    // Add to the master list
     237                    pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
     238                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     239                    xMaster->data.F32[numMaster] = xImage->data.F32[j];
     240                    yMaster->data.F32[numMaster] = yImage->data.F32[j];
     241                    matches->data[numMaster] = match;
     242                    numMaster++;
     243                    xMaster->n = yMaster->n = matches->n = numMaster;
     244                }
     245
     246            }
     247            psFree(coords);
     248            psFree(tree);
     249        }
     250
    258251        psFree(boundsImage);
    259252        psFree(xImage);
     
    263256    }
    264257
    265     // Now cull the matches that contain only a single star
    266     int numGood = 0;                    // Number of good matches
     258    if (cullSingles) {
     259        // Now cull the matches that contain only a single star
     260        int numGood = 0;                    // Number of good matches
     261        for (int i = 0; i < matches->n; i++) {
     262            pmSourceMatch *match = matches->data[i]; // Match of interest
     263            if (match->num > 1) {
     264                if (i != numGood) {
     265                    psFree(matches->data[numGood]);
     266                    matches->data[numGood] = psMemIncrRefCounter(match);
     267                }
     268                numGood++;
     269            }
     270        }
     271        matches->n = numGood;
     272        for (int i = numGood; i < numMaster; i++) {
     273            psFree(matches->data[i]);
     274            matches->data[i] = NULL;
     275        }
     276    }
     277
     278    return matches;
     279}
     280
     281
     282psArray *pmSourceMatchMerge(psArray *sourceArrays, float radius)
     283{
     284    PS_ASSERT_ARRAY_NON_NULL(sourceArrays, NULL);
     285    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, NULL);
     286
     287    psArray *matches = pmSourceMatchSources(sourceArrays, radius, false); // Source matches
     288    if (!matches) {
     289        psError(PS_ERR_UNKNOWN, false, "Unable to match source lists.");
     290        return NULL;
     291    }
     292
     293    int index = 0;                      // Index to current position on list
    267294    for (int i = 0; i < matches->n; i++) {
     295        pmSource *source = NULL;        // Source to put in merged list
    268296        pmSourceMatch *match = matches->data[i]; // Match of interest
    269         if (match->num > 1) {
    270             if (i != numGood) {
    271                 psFree(matches->data[numGood]);
    272                 matches->data[numGood] = psMemIncrRefCounter(match);
    273                 //                psFree(match);
    274             }
    275             numGood++;
    276         }
    277     }
    278     matches->n = numGood;
    279     for (int i = numGood; i < numMaster; i++) {
    280         psFree(matches->data[i]);
    281         matches->data[i] = NULL;
    282     }
     297        for (int j = 0; j < match->num && !source; j++) {
     298            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) {
     299                continue;
     300            }
     301            int imgIndex = match->image->data.S32[j]; // Index of image
     302            int srcIndex = match->index->data.S32[j]; // Index of source for image
     303            psArray *list = sourceArrays->data[imgIndex]; // List of interest
     304            source = list->data[srcIndex];
     305            break;
     306        }
     307
     308        if (source) {
     309            psFree(matches->data[index]);
     310            matches->data[index] = psMemIncrRefCounter(source);
     311            index++;
     312        }
     313    }
     314
     315    // Clear out the rest of the list
     316    int num = index;                    // Number of good sources
     317    for (; index < matches->n; index++) {
     318        psFree(matches->data[index]);
     319        matches->data[index] = NULL;
     320    }
     321    matches->n = num;
    283322
    284323    return matches;
    285324}
    286 
    287325
    288326// Iterate on the star magnitudes and image transparencies
     
    290328static float sourceMatchRelphotIterate(psVector *trans, // Transparencies
    291329                                       psVector *stars, // Star magnitudes
     330                                       psVector *badImage, // Bad image mask
    292331                                       const psArray *matches, // Array of matches
    293332                                       const psVector *zp, // Zero points for each image (incl. airmass term)
     
    363402    for (int i = 0; i < numImages; i++) {
    364403        trans->data.F32[i] = accum->data.F64[i] / accumErr->data.F64[i];
    365 
     404        if (!isfinite(trans->data.F32[i])) {
     405            badImage->data.U8[i] = 0xFF;
     406        }
    366407        psTrace("psModules.objects", 3, "Transparency for image %d: %f\n", i, trans->data.F32[i]);
    367408    }
     
    379420            }
    380421            int index = match->image->data.U32[j]; // Image index
     422            if (badImage->data.U8[index]) {
     423                continue;
     424            }
    381425            float mag = match->mag->data.F32[j]; // Measured magnitude
    382426            float magErr2 = PS_SQR(match->magErr->data.F32[j]) + sysErr2; // Error in measured magnitude
     
    402446static int sourceMatchRelphotPhotometric(psVector *photo, // Photometric determination
    403447                                         const psVector *trans, // Estimated transparencies
     448                                         const psVector *badImage, // Bad image?
    404449                                         int transIter, // Iterations for transparency
    405450                                         float transClip, // Clipping level for transparency
     
    409454    psAssert(photo && photo->type.type == PS_TYPE_U8, "Need photometric determination");
    410455    psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies");
     456    psAssert(badImage && badImage->type.type == PS_TYPE_U8, "Need bad image determination");
    411457
    412458    int numImages = photo->n;              // Number of images
    413459
    414460    psAssert(trans->n == numImages, "Not enough transparencies: %ld", trans->n);
     461    psAssert(badImage->n == numImages, "Not enough bad image determinations: %ld", badImage->n);
    415462    psAssert(transIter >= 0, "Iterations for transparency must be non-negative: %d", transIter);
    416463    psAssert(transClip > 0, "Clipping level for transparency must be positive: %f", transClip);
     
    421468    stats->clipSigma = transClip;
    422469
    423     if (!psVectorStats(stats, trans, NULL, NULL, 0)) {
     470    if (!psVectorStats(stats, trans, NULL, badImage, 0xFF)) {
    424471        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies.");
    425472        psFree(stats);
     
    432479    int numPhoto = 0;                   // Number of photometric images
    433480    for (int i = 0; i < numImages; i++) {
     481        if (badImage->data.U8[i]) {
     482            continue;
     483        }
    434484        if (trans->data.F32[i] < thresh) {
    435485            photo->data.U8[i] = 0xFF;
     
    451501                                      const psVector *zp, // Zero points for each image
    452502                                      const psVector *photo, // Photometric image?
     503                                      const psVector *badImage, // Bad image?
    453504                                      float starClip, // Clipping for stars
    454505                                      float sysErr2 // Systematic error squared
     
    468519    psAssert(!photo || photo->type.type == PS_TYPE_U8, "Photometric determination is wrong type");
    469520    psAssert(!photo || photo->n == numImages, "Not enough photometric determinations: %ld", photo->n);
     521    psAssert(!badImage || badImage->type.type == PS_TYPE_U8, "Photometric determination is wrong type");
     522    psAssert(!badImage || badImage->n == numImages, "Not enough bad determinations: %ld", badImage->n);
    470523
    471524    starClip = PS_SQR(starClip);
     
    481534            numMeasurements++;
    482535            int index = match->image->data.U32[j]; // Image index
     536            if (badImage->data.U8[index]) {
     537                continue;
     538            }
    483539            float mag = match->mag->data.F32[j]; // Measured magnitude
    484540            float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
     
    489545            float dev = mag + cal - stars->data.F32[i]; // Deviation
    490546
    491             // only reject detections from photometric images (non-photometric images can
    492             // have large errors.  XXX Or: allow a much higher rejection threshold
    493             if (photo->data.U8[index]) {
    494               if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
    495                 numRejected++;
    496                 match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
    497               }
    498             }
     547            // only reject detections from photometric images (non-photometric images can
     548            // have large errors.  XXX Or: allow a much higher rejection threshold
     549            if (photo->data.U8[index]) {
     550                if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
     551                    numRejected++;
     552                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
     553                }
     554            }
    499555        }
    500556    }
     
    506562psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches
    507563                               const psVector *zp, // Zero points for each image (including airmass term)
    508                                int maxIter, // Maximum number of iterations
    509564                               float tol, // Relative tolerance for convergence
     565                               int iter1, // Number of iterations for pass 1
     566                               float rej1, // Limit on rejection between iterations for pass 1
     567                               float sys1, // Systematic error in measurements for pass 1
     568                               int iter2, // Number of iterations for pass 2
     569                               float rej2, // Limit on rejection between iterations for pass 2
     570                               float sys2, // Systematic error in measurements for pass 2
    510571                               float rejLimit, // Limit on rejection between iterations
    511572                               int transIter, // Clipping iterations for transparency
    512573                               float transClip, // Clipping level for transparency
    513                                float photoLevel, // Level at which we declare image is photometric
    514                                float starClip, // Clipping for stars
    515                                float sysErr // Systematic error in measurements
     574                               float photoLevel // Level at which we declare image is photometric
    516575                               )
    517576{
     
    521580    PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL);
    522581
    523     sysErr *= sysErr;
     582    sys1 *= sys1;
     583    sys2 *= sys2;
    524584
    525585    int numImages = zp->n;              // Number of images
     
    529589    psVector *photo = psVectorAlloc(numImages, PS_TYPE_U8); // Photometric determination for each image
    530590    psVectorInit(photo, 0);
     591    psVector *badImage = psVectorAlloc(numImages, PS_TYPE_U8); // Bad image?
     592    psVectorInit(badImage, 0);
    531593    psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
    532594
    533     float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, sysErr); // chi^2 for solution
     595    float chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp,
     596                                           photo, sys1); // chi^2 for solution
    534597    psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
    535598    float lastChi2 = INFINITY;          // chi^2 on last iteration
    536599    float fracRej = INFINITY;        // Fraction of measurements rejected
    537600
    538     // in the first passes, the transparencies are not well deteremined: use high systematic error and the rejection thresholds
    539     for (int i = 0; i < 5; i++) {
     601    // In the first passes, the transparencies are not well deteremined: use high systematic error and
     602    // rejection thresholds
     603    for (int i = 0; i < iter1; i++) {
    540604
    541605        // Identify photometric nights
    542         int numPhoto = sourceMatchRelphotPhotometric(photo, trans, transIter, transClip, photoLevel); // Number of photometric images
    543         if (numPhoto < 0) {
    544             psError(PS_ERR_UNKNOWN, false, "Unable to perform photometric determination");
    545             psFree(trans);
    546             psFree(photo);
    547             psFree(stars);
    548             return NULL;
    549         }
    550         psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
    551 
    552         // XXX use 20 sigma rejection and 0.1 mag systematic error (move these to the recipe)
    553         fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, 20.0, PS_SQR(0.1));
    554         psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
    555 
    556         // XXX use 0.05 mag systematic error (move these to the recipe)
    557         chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, PS_SQR(0.1)); // chi^2 for solution
    558         psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
    559     }
    560 
    561     for (int i = 0; i < maxIter && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
    562         lastChi2 = chi2;
    563 
    564         // Identify photometric nights
    565         int numPhoto = sourceMatchRelphotPhotometric(photo, trans, transIter, transClip,
     606        int numPhoto = sourceMatchRelphotPhotometric(photo, trans, badImage, transIter, transClip,
    566607                                                     photoLevel); // Number of photometric images
    567608        if (numPhoto < 0) {
     
    572613            return NULL;
    573614        }
    574         psTrace("psModules.objects", 3, "Determined %d/%d are photometric", numPhoto, numImages);
    575 
    576         fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, starClip, sysErr);
    577         psTrace("psModules.objects", 3, "%f%% of measurements rejected", fracRej * 100);
    578 
    579         chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, sysErr); // chi^2 for solution
    580         psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     615        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
     616
     617        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1);
     618        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
     619
     620        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1);
     621        psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     622    }
     623
     624    for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
     625        lastChi2 = chi2;
     626
     627        // Identify photometric nights
     628        int numPhoto = sourceMatchRelphotPhotometric(photo, trans, badImage, transIter, transClip,
     629                                                     photoLevel); // Number of photometric images
     630        if (numPhoto < 0) {
     631            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometric determination");
     632            psFree(trans);
     633            psFree(photo);
     634            psFree(stars);
     635            return NULL;
     636        }
     637        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages);
     638
     639        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2);
     640        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100);
     641
     642        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2);
     643        psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
    581644    }
    582645
    583646    psFree(photo);
     647    psFree(badImage);
    584648    psFree(stars);
    585649
Note: See TracChangeset for help on using the changeset viewer.