Changeset 26532
- Timestamp:
- Jan 7, 2010, 12:05:55 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/20091201/psModules/src/astrom/pmAstrometryObjects.c
r26260 r26532 108 108 } 109 109 110 if (found1->data.S8[i]) {111 i++;112 continue;113 }114 if (found2->data.S8[j]) {115 j++;116 continue;117 }110 if (found1->data.S8[i]) { 111 i++; 112 continue; 113 } 114 if (found2->data.S8[j]) { 115 j++; 116 continue; 117 } 118 118 119 119 jStart = j; … … 128 128 continue; 129 129 } 130 if (found2->data.S8[j]) {131 j++;132 continue;133 }130 if (found2->data.S8[j]) { 131 j++; 132 continue; 133 } 134 134 135 135 // got a match; add to output list … … 138 138 psFree (match); 139 139 140 found1->data.S8[i] = 1;141 found2->data.S8[j] = 1;140 found1->data.S8[i] = 1; 141 found2->data.S8[j] = 1; 142 142 143 143 j++; … … 196 196 psArray *matches = match_lists(x1, y1, x2, y2, sorted1, sorted2, RADIUS); \ 197 197 \ 198 psFree( (void *)sorted1); \199 psFree( (void *)sorted2); \198 psFree(sorted1); \ 199 psFree(sorted2); \ 200 200 psFree(x1); \ 201 201 psFree(y1); \ … … 314 314 if (!xRes) abort(); 315 315 psFree (xFit); 316 316 317 317 psVector *yFit = psPolynomial2DEvalVector (map->y, X, Y); 318 318 if (!yFit) abort(); … … 322 322 323 323 // extract a high-quality subset (unmasked, S/N > XXX) and position errors 324 // XXX for now, generate a position error based on the magnitude error 324 // XXX for now, generate a position error based on the magnitude error 325 325 psVector *xErr = psVectorAllocEmpty (match->n, PS_TYPE_F32); 326 326 psVector *yErr = psVectorAllocEmpty (match->n, PS_TYPE_F32); … … 329 329 330 330 for (int i = 0; i < match->n; i++) { 331 if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;331 if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; 332 332 pmAstromMatch *pair = match->data[i]; 333 333 pmAstromObj *rawStar = raw->data[pair->raw]; … … 335 335 if (rawStar->dMag > 0.02) continue; 336 336 337 // two likely failure values: NAN or 0.0 --> use dMag in this case338 float xErrValue, yErrValue;339 if (isfinite(rawStar->chip->xErr) && (rawStar->chip->xErr > 0.0)) {340 xErrValue = rawStar->chip->xErr;341 } else {342 xErrValue = PS_MAX(0.005, rawStar->dMag);343 }344 if (isfinite(rawStar->chip->yErr) && (rawStar->chip->yErr > 0.0)) {345 yErrValue = rawStar->chip->yErr;346 } else {347 yErrValue = PS_MAX(0.005, rawStar->dMag);348 }349 350 psVectorAppend (xErr, xErrValue);351 psVectorAppend (yErr, yErrValue);352 psVectorAppend (xResGood, xRes->data.F32[i]);353 psVectorAppend (yResGood, yRes->data.F32[i]);337 // two likely failure values: NAN or 0.0 --> use dMag in this case 338 float xErrValue, yErrValue; 339 if (isfinite(rawStar->chip->xErr) && (rawStar->chip->xErr > 0.0)) { 340 xErrValue = rawStar->chip->xErr; 341 } else { 342 xErrValue = PS_MAX(0.005, rawStar->dMag); 343 } 344 if (isfinite(rawStar->chip->yErr) && (rawStar->chip->yErr > 0.0)) { 345 yErrValue = rawStar->chip->yErr; 346 } else { 347 yErrValue = PS_MAX(0.005, rawStar->dMag); 348 } 349 350 psVectorAppend (xErr, xErrValue); 351 psVectorAppend (yErr, yErrValue); 352 psVectorAppend (xResGood, xRes->data.F32[i]); 353 psVectorAppend (yResGood, yRes->data.F32[i]); 354 354 } 355 355 … … 434 434 // Get the minimum and maximum values 435 435 if (!psVectorStats(stats, myVector, NULL, NULL, 0)) { 436 psFree(stats);437 return NAN;436 psFree(stats); 437 return NAN; 438 438 } 439 439 min = stats->min; 440 440 max = stats->max; 441 441 if (isnan(min) || isnan(max)) { 442 psFree(stats);443 return NAN;442 psFree(stats); 443 return NAN; 444 444 } 445 445 … … 448 448 // If all data points have the same value, then we set the appropriate members of stats and return. 449 449 if (fabs(max - min) <= FLT_EPSILON) { 450 psFree (stats);451 return 0.0;452 } 453 454 // Define the histogram bin size. 450 psFree (stats); 451 return 0.0; 452 } 453 454 // Define the histogram bin size. 455 455 float binSize = 0.001; 456 456 long numBins = PS_MIN(100000, (max - min) / binSize); // Number of bins … … 463 463 464 464 if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) { 465 // if psVectorHistogram returns false, we have a programming error466 psAbort ("Unable to generate histogram");465 // if psVectorHistogram returns false, we have a programming error 466 psAbort ("Unable to generate histogram"); 467 467 } 468 468 if (psTraceGetLevel("psModules.astrom") >= 8) { 469 PS_VECTOR_PRINT_F32(histogram->bounds);470 PS_VECTOR_PRINT_F32(histogram->nums);469 PS_VECTOR_PRINT_F32(histogram->bounds); 470 PS_VECTOR_PRINT_F32(histogram->nums); 471 471 } 472 472 … … 475 475 cumulative->nums->data.F32[0] = histogram->nums->data.F32[0]; 476 476 for (long i = 1; i < histogram->nums->n; i++) { 477 cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];478 cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];477 cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i]; 478 cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i]; 479 479 } 480 480 if (psTraceGetLevel("psModules.astrom") >= 8) { 481 PS_VECTOR_PRINT_F32(cumulative->bounds);482 PS_VECTOR_PRINT_F32(cumulative->nums);481 PS_VECTOR_PRINT_F32(cumulative->bounds); 482 PS_VECTOR_PRINT_F32(cumulative->nums); 483 483 } 484 484 … … 837 837 } 838 838 839 if (psTraceGetLevel("psModules.astrom") >= 5) { 840 char line[16]; 841 psFits *fits = psFitsOpen ("grid.image.fits", "w"); 842 psFitsWriteImage (fits, NULL, gridNP, 0, NULL); 843 psFitsClose (fits); 844 fprintf (stderr, "wrote grid image, press return to continue\n"); 845 fgets (line, 15, stdin); 846 } 839 if (psTraceGetLevel("psModules.astrom") >= 5) { 840 char line[16]; 841 psFits *fits = psFitsOpen ("grid.image.fits", "w"); 842 psFitsWriteImage (fits, NULL, gridNP, 0, NULL); 843 psFitsClose (fits); 844 fprintf (stderr, "wrote grid image, press return to continue\n"); 845 if (!fgets (line, 15, stdin)) { 846 fprintf(stderr, "Error waiting for RETURN."); 847 } 848 } 847 849 848 850 // only check bins with at least 1/2 of max bin … … 1189 1191 psArray *pmAstromRadiusMatchUniq (psArray *rawstars, psArray *refstars, psArray *matches) { 1190 1192 1191 // I have the matches between the refstars and the rawstars. 1193 // I have the matches between the refstars and the rawstars. 1192 1194 // For each refstar, find the single match which has the smallest radius and reject 1193 // all others. 1195 // all others. 1194 1196 1195 1197 // create an array of refstars->n arrays, each containing all of the matches for the 1196 // given refstar. 1198 // given refstar. 1197 1199 1198 1200 psArray *refstarMatches = psArrayAlloc (refstars->n); … … 1200 1202 for (int i = 0; i < matches->n; i++) { 1201 1203 1202 pmAstromMatch *match = matches->data[i];1203 1204 // accumulate this refstar match on the array for this refstar (create if needed)1205 psArray *refSet = refstarMatches->data[match->ref];1206 if (!refSet) {1207 refstarMatches->data[match->ref] = psArrayAllocEmpty (8);1208 refSet = refstarMatches->data[match->ref];1209 }1210 1211 pmAstromMatchInfo *matchInfo = pmAstromMatchInfoAlloc();1212 1213 pmAstromObj *refStar = refstars->data[match->ref];1214 pmAstromObj *rawStar = rawstars->data[match->raw];1215 1216 matchInfo->match = match; // reference to the match of interest1217 matchInfo->radius = hypot (refStar->FP->x - rawStar->FP->x, refStar->FP->y - rawStar->FP->y);1218 1219 psArrayAdd (refSet, 8, matchInfo); // matchInfo->match is just a reference1220 psFree (matchInfo);1204 pmAstromMatch *match = matches->data[i]; 1205 1206 // accumulate this refstar match on the array for this refstar (create if needed) 1207 psArray *refSet = refstarMatches->data[match->ref]; 1208 if (!refSet) { 1209 refstarMatches->data[match->ref] = psArrayAllocEmpty (8); 1210 refSet = refstarMatches->data[match->ref]; 1211 } 1212 1213 pmAstromMatchInfo *matchInfo = pmAstromMatchInfoAlloc(); 1214 1215 pmAstromObj *refStar = refstars->data[match->ref]; 1216 pmAstromObj *rawStar = rawstars->data[match->raw]; 1217 1218 matchInfo->match = match; // reference to the match of interest 1219 matchInfo->radius = hypot (refStar->FP->x - rawStar->FP->x, refStar->FP->y - rawStar->FP->y); 1220 1221 psArrayAdd (refSet, 8, matchInfo); // matchInfo->match is just a reference 1222 psFree (matchInfo); 1221 1223 } 1222 1224 … … 1227 1229 for (int i = 0; i < refstars->n; i++) { 1228 1230 1229 psArray *refSet = refstarMatches->data[i];1230 if (!refSet) continue;1231 if (refSet->n == 0) continue; // not certain how this can happen...1232 1233 if (refSet->n == 1) {1234 pmAstromMatchInfo *matchInfo = refSet->data[0];1235 psArrayAdd (unique, 32, matchInfo->match);1236 continue;1237 }1238 1239 pmAstromMatchInfo *matchInfo = refSet->data[0];1240 float minRadius = matchInfo->radius;1241 pmAstromMatch *minMatch = matchInfo->match;1242 for (int j = 1; j < refSet->n; j++) {1243 pmAstromMatchInfo *matchInfo = refSet->data[j];1244 if (minRadius < matchInfo->radius) continue;1245 minMatch = matchInfo->match;1246 minRadius = matchInfo->radius;1247 }1248 1249 psArrayAdd (unique, 32, minMatch); // minMatch is just a reference to a match on matches,1231 psArray *refSet = refstarMatches->data[i]; 1232 if (!refSet) continue; 1233 if (refSet->n == 0) continue; // not certain how this can happen... 1234 1235 if (refSet->n == 1) { 1236 pmAstromMatchInfo *matchInfo = refSet->data[0]; 1237 psArrayAdd (unique, 32, matchInfo->match); 1238 continue; 1239 } 1240 1241 pmAstromMatchInfo *matchInfo = refSet->data[0]; 1242 float minRadius = matchInfo->radius; 1243 pmAstromMatch *minMatch = matchInfo->match; 1244 for (int j = 1; j < refSet->n; j++) { 1245 pmAstromMatchInfo *matchInfo = refSet->data[j]; 1246 if (minRadius < matchInfo->radius) continue; 1247 minMatch = matchInfo->match; 1248 minRadius = matchInfo->radius; 1249 } 1250 1251 psArrayAdd (unique, 32, minMatch); // minMatch is just a reference to a match on matches, 1250 1252 } 1251 1253
Note:
See TracChangeset
for help on using the changeset viewer.
