- Timestamp:
- Mar 17, 2009, 12:07:42 PM (17 years ago)
- Location:
- branches/cnb_branches/cnb_branch_20090301/psModules
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/objects/pmSourceMatch.c (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/cnb_branches/cnb_branch_20090301/psModules
-
Property svn:mergeinfo
set to (toggle deleted branches)
/branches/cnb_branches/cnb_branch_20090215/psModules merged eligible /trunk/psModules merged eligible /branches/cnb_branch_20090215/psModules 21495-22685 /branches/eam_branches/eam_branch_20090303/psModules 23158-23228
-
Property svn:mergeinfo
set to (toggle deleted branches)
-
branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmSourceMatch.c
r21266 r23351 16 16 #define SOURCE_FAINTEST 50.0 // Faintest magnitude to consider 17 17 #define SOURCES_MAX_LEAF 2 // Maximum number of points on a tree leaf 18 18 #define ARRAY_BUFFER 16 // Buffer for array 19 19 20 20 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 63 63 pmSource *source = sources->data[i]; // Source of interest 64 64 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; 69 69 70 70 float xSrc, ySrc; // Coordinates of source … … 113 113 } 114 114 115 pmSourceMatch *pmSourceMatchAlloc(int num // Maximum number of images 116 ) 115 pmSourceMatch *pmSourceMatchAlloc(void) 117 116 { 118 117 pmSourceMatch *match = psAlloc(sizeof(pmSourceMatch)); // Match data … … 120 119 121 120 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); 127 126 128 127 return match; … … 136 135 { 137 136 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); 142 143 143 144 match->mag->data.F32[num] = mag; … … 148 149 match->num++; 149 150 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 155 psArray *pmSourceMatchSources(const psArray *sourceArrays, float radius, bool cullSingles) 155 156 { 156 157 PS_ASSERT_ARRAY_NON_NULL(sourceArrays, NULL); … … 185 186 matches = psArrayAlloc(numSources); 186 187 for (int j = 0; j < numSources; j++) { 187 pmSourceMatch *match = pmSourceMatchAlloc( numImages); // Match data188 pmSourceMatch *match = pmSourceMatchAlloc(); // Match data 188 189 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); 189 190 matches->data[j] = match; 190 191 } 191 192 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 197 196 psTrace("psModules.objects", 7, "Bounds don't overlap\n"); 198 // Add everything in to the master list199 197 long size = numMaster + numSources; // New size 200 198 xMaster = psVectorRealloc(xMaster, size); … … 207 205 numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); 208 206 for (int j = 0, k = numMaster; j < numSources; j++, k++) { 209 pmSourceMatch *match = pmSourceMatchAlloc( numImages); // Match data207 pmSourceMatch *match = pmSourceMatchAlloc(); // Match data 210 208 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); 211 209 matches->data[k] = match; … … 215 213 yMaster->n = size; 216 214 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 258 251 psFree(boundsImage); 259 252 psFree(xImage); … … 263 256 } 264 257 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 282 psArray *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 267 294 for (int i = 0; i < matches->n; i++) { 295 pmSource *source = NULL; // Source to put in merged list 268 296 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; 283 322 284 323 return matches; 285 324 } 286 287 325 288 326 // Iterate on the star magnitudes and image transparencies … … 290 328 static float sourceMatchRelphotIterate(psVector *trans, // Transparencies 291 329 psVector *stars, // Star magnitudes 330 psVector *badImage, // Bad image mask 292 331 const psArray *matches, // Array of matches 293 332 const psVector *zp, // Zero points for each image (incl. airmass term) … … 363 402 for (int i = 0; i < numImages; i++) { 364 403 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 } 366 407 psTrace("psModules.objects", 3, "Transparency for image %d: %f\n", i, trans->data.F32[i]); 367 408 } … … 379 420 } 380 421 int index = match->image->data.U32[j]; // Image index 422 if (badImage->data.U8[index]) { 423 continue; 424 } 381 425 float mag = match->mag->data.F32[j]; // Measured magnitude 382 426 float magErr2 = PS_SQR(match->magErr->data.F32[j]) + sysErr2; // Error in measured magnitude … … 402 446 static int sourceMatchRelphotPhotometric(psVector *photo, // Photometric determination 403 447 const psVector *trans, // Estimated transparencies 448 const psVector *badImage, // Bad image? 404 449 int transIter, // Iterations for transparency 405 450 float transClip, // Clipping level for transparency … … 409 454 psAssert(photo && photo->type.type == PS_TYPE_U8, "Need photometric determination"); 410 455 psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies"); 456 psAssert(badImage && badImage->type.type == PS_TYPE_U8, "Need bad image determination"); 411 457 412 458 int numImages = photo->n; // Number of images 413 459 414 460 psAssert(trans->n == numImages, "Not enough transparencies: %ld", trans->n); 461 psAssert(badImage->n == numImages, "Not enough bad image determinations: %ld", badImage->n); 415 462 psAssert(transIter >= 0, "Iterations for transparency must be non-negative: %d", transIter); 416 463 psAssert(transClip > 0, "Clipping level for transparency must be positive: %f", transClip); … … 421 468 stats->clipSigma = transClip; 422 469 423 if (!psVectorStats(stats, trans, NULL, NULL, 0)) {470 if (!psVectorStats(stats, trans, NULL, badImage, 0xFF)) { 424 471 psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies."); 425 472 psFree(stats); … … 432 479 int numPhoto = 0; // Number of photometric images 433 480 for (int i = 0; i < numImages; i++) { 481 if (badImage->data.U8[i]) { 482 continue; 483 } 434 484 if (trans->data.F32[i] < thresh) { 435 485 photo->data.U8[i] = 0xFF; … … 451 501 const psVector *zp, // Zero points for each image 452 502 const psVector *photo, // Photometric image? 503 const psVector *badImage, // Bad image? 453 504 float starClip, // Clipping for stars 454 505 float sysErr2 // Systematic error squared … … 468 519 psAssert(!photo || photo->type.type == PS_TYPE_U8, "Photometric determination is wrong type"); 469 520 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); 470 523 471 524 starClip = PS_SQR(starClip); … … 481 534 numMeasurements++; 482 535 int index = match->image->data.U32[j]; // Image index 536 if (badImage->data.U8[index]) { 537 continue; 538 } 483 539 float mag = match->mag->data.F32[j]; // Measured magnitude 484 540 float magErr = match->magErr->data.F32[j]; // Error in measured magnitude … … 489 545 float dev = mag + cal - stars->data.F32[i]; // Deviation 490 546 491 // only reject detections from photometric images (non-photometric images can492 // have large errors. XXX Or: allow a much higher rejection threshold493 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 } 499 555 } 500 556 } … … 506 562 psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches 507 563 const psVector *zp, // Zero points for each image (including airmass term) 508 int maxIter, // Maximum number of iterations509 564 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 510 571 float rejLimit, // Limit on rejection between iterations 511 572 int transIter, // Clipping iterations for transparency 512 573 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 516 575 ) 517 576 { … … 521 580 PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL); 522 581 523 sysErr *= sysErr; 582 sys1 *= sys1; 583 sys2 *= sys2; 524 584 525 585 int numImages = zp->n; // Number of images … … 529 589 psVector *photo = psVectorAlloc(numImages, PS_TYPE_U8); // Photometric determination for each image 530 590 psVectorInit(photo, 0); 591 psVector *badImage = psVectorAlloc(numImages, PS_TYPE_U8); // Bad image? 592 psVectorInit(badImage, 0); 531 593 psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star 532 594 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 534 597 psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2); 535 598 float lastChi2 = INFINITY; // chi^2 on last iteration 536 599 float fracRej = INFINITY; // Fraction of measurements rejected 537 600 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++) { 540 604 541 605 // 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, 566 607 photoLevel); // Number of photometric images 567 608 if (numPhoto < 0) { … … 572 613 return NULL; 573 614 } 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); 581 644 } 582 645 583 646 psFree(photo); 647 psFree(badImage); 584 648 psFree(stars); 585 649
Note:
See TracChangeset
for help on using the changeset viewer.
