- Timestamp:
- Jan 28, 2009, 1:36:37 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/cnb_branch_20090113/psastro/src/psastroAstromGuess.c
r20805 r21208 29 29 psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE); 30 30 if (!recipe) { 31 32 31 psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); 32 return false; 33 33 } 34 34 … … 36 36 bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); 37 37 if (!status) { 38 38 useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); 39 39 } 40 40 … … 42 42 pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); 43 43 if (!input) { 44 45 44 psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); 45 return false; 46 46 } 47 47 … … 49 49 double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); 50 50 if (!status) { 51 psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); 52 return false; 53 } 51 psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); 52 return false; 53 } 54 54 55 55 psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32); … … 67 67 bool bilevelAstrometry = false; 68 68 if (!useModel) { 69 69 psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); 70 70 } 71 71 … … 75 75 if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } 76 76 77 78 79 77 if (!useModel) { 78 if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; 79 } 80 80 81 81 if (newFPA) { 82 82 newFPA = false; 83 84 83 while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; 84 while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; 85 85 RAminSky = fpa->toSky->R - M_PI; 86 86 RAmaxSky = fpa->toSky->R + M_PI; 87 87 } 88 88 89 90 { 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 89 // report and save the current best guess for the chip 0,0 pixel coordinates 90 { 91 psPlane ptCH, ptFP, ptTP; 92 psSphere ptSky; 93 94 ptCH.x = 0; 95 ptCH.y = 0; 96 psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); 97 psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); 98 psDeproject (&ptSky, &ptTP, fpa->toSky); 99 psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); 100 101 psVectorAppend (cornerL, ptFP.x); 102 psVectorAppend (cornerM, ptFP.y); 103 psVectorAppend (cornerP, ptTP.x); 104 psVectorAppend (cornerQ, ptTP.y); 105 psVectorAppend (cornerR, ptSky.r); 106 psVectorAppend (cornerD, ptSky.d); 107 } 108 108 109 109 // apply the new WCS guess data to all of the data in the readouts … … 119 119 if (rawstars == NULL) { continue; } 120 120 121 121 *nStars += rawstars->n; 122 122 for (int i = 0; i < rawstars->n; i++) { 123 123 pmAstromObj *raw = rawstars->data[i]; … … 138 138 } 139 139 140 141 142 143 144 145 p sastroVisualPlotRawStars(rawstars, fpa, chip, recipe);146 147 148 149 140 // dump or plot the resulting projected positions 141 if (psTraceGetLevel("psastro.dump") > 0) { 142 psastroDumpRawstars (rawstars, fpa, chip); 143 } 144 145 pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe); 146 147 if (psTraceGetLevel("psastro.plot") > 0) { 148 psastroPlotRawstars (rawstars, fpa, chip, recipe); 149 } 150 150 } 151 151 } … … 157 157 psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR", PS_META_REPLACE, "", *nStars); 158 158 if (*nStars == 0) { 159 160 161 162 } 163 164 psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 165 DEG_RAD*RAmin, DEG_RAD*DECmin, 166 159 psLogMsg ("psastro", 2, "no sources available for astrometry\n"); 160 psFree (view); 161 return true; 162 } 163 164 psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 165 DEG_RAD*RAmin, DEG_RAD*DECmin, 166 DEG_RAD*RAmax, DEG_RAD*DECmax); 167 167 168 168 psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", RAmin); … … 203 203 pmHDU *hdu = pmFPAviewThisHDU (view, fpa); 204 204 if (bilevelAstrometry) { 205 206 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 207 208 } 205 if (!pmAstromReadBilevelChip (chip, hdu->header)) { 206 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 207 return false; 208 } 209 209 } else { 210 211 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 212 213 } 210 if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { 211 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 212 return false; 213 } 214 214 } 215 215 return true; … … 225 225 // load mosaic-level astrometry? 226 226 if (phu) { 227 228 229 230 227 char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); 228 if (ctype) { 229 *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); 230 } 231 231 } 232 232 if (*bilevelAstrometry) { 233 234 } 233 pmAstromReadBilevelMosaic (fpa, phu->header); 234 } 235 235 psFree (view); 236 236 return true; … … 245 245 pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); 246 246 if (!input) { 247 248 247 psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); 248 return false; 249 249 } 250 250 … … 277 277 if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } 278 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 279 // XXX we are currently inconsistent with marking the good vs the bad data 280 // psastroChipAstrom sets data_exists to false if the fit is bad. this is 281 // probably wrong since it implies there is no data! 282 283 // skip chips for which the astrometry failed (NASTRO == 0) 284 if (!chip->cells->n) goto skip_chip; 285 pmCell *cell = chip->cells->data[0]; 286 if (!cell) goto skip_chip; 287 288 if (!cell->readouts->n) goto skip_chip; 289 pmReadout *readout = cell->readouts->data[0]; 290 if (!readout) goto skip_chip; 291 292 psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); 293 if (!updates) goto skip_chip; 294 295 int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); 296 if (!nAstro) goto skip_chip; 297 298 float astError = psMetadataLookupF32 (&status, updates, "CERROR"); 299 if (fabs(astError) < 1e-6) goto skip_chip; 300 301 psPlane ptCH, ptFP, ptTP; 302 psSphere ptSky; 303 304 ptCH.x = 0; 305 ptCH.y = 0; 306 psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); 307 psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); 308 psDeproject (&ptSky, &ptTP, fpa->toSky); 309 psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); 310 311 // new corner locations based on the calibrated astrometry 312 psVectorAppend (cornerLn, ptFP.x); 313 psVectorAppend (cornerMn, ptFP.y); 314 psVectorAppend (cornerPn, ptTP.x); 315 psVectorAppend (cornerQn, ptTP.y); 316 psVectorAppend (cornerRn, ptSky.r); 317 psVectorAppend (cornerDn, ptSky.d); 318 psVectorAppend (cornerMK, 0); 319 continue; 320 320 321 321 skip_chip: 322 323 324 325 326 327 328 329 322 // new corner locations based on the calibrated astrometry 323 psVectorAppend (cornerLn, 0.0); 324 psVectorAppend (cornerMn, 0.0); 325 psVectorAppend (cornerPn, 0.0); 326 psVectorAppend (cornerQn, 0.0); 327 psVectorAppend (cornerRn, 0.0); 328 psVectorAppend (cornerDn, 0.0); 329 psVectorAppend (cornerMK, 1); 330 330 } 331 331 … … 336 336 337 337 for (int i = 0; i < cornerRo->n; i++) { 338 339 340 341 342 343 344 345 346 347 338 339 psPlane ptTP; 340 psSphere ptSky; 341 342 ptSky.r = cornerRo->data.F32[i]; 343 ptSky.d = cornerDo->data.F32[i]; 344 345 psProject (&ptTP, &ptSky, fpa->toSky); 346 psVectorAppend (cornerPs, ptTP.x); 347 psVectorAppend (cornerQs, ptTP.y); 348 348 } 349 349 … … 351 351 map->x->coeffMask[1][1] = PS_POLY_MASK_SET; 352 352 map->y->coeffMask[1][1] = PS_POLY_MASK_SET; 353 353 354 354 // fit the valid chips, mask the invalid chips 355 355 psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs); 356 356 psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs); 357 357 358 358 // apply the linear fit... 359 359 psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs); … … 364 364 psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf); 365 365 366 p sastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);366 pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); 367 367 368 368 psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); … … 374 374 float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]); 375 375 float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]); 376 376 377 377 psLogMsg ("psastro", 3, "boresite offset : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]); 378 378 psLogMsg ("psastro", 3, "boresite angle : %f, scale: %f", angle*PS_DEG_RAD, scale); … … 382 382 psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER"); 383 383 if (!header) { 384 385 386 384 header = psMetadataAlloc(); 385 psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); 386 psFree (header); // drop this reference 387 387 } 388 388 … … 395 395 396 396 if (DEBUG) { 397 398 399 400 cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], 401 402 403 397 FILE *f = fopen ("corners.dat", "w"); 398 for (int i = 0; i < cornerRo->n; i++) { 399 fprintf (f, "%10.6f %10.6f %9.2f %9.2f %9.2f %9.2f | %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f\n", 400 cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], 401 cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); 402 } 403 fclose (f); 404 404 } 405 405 … … 424 424 psFree (map); 425 425 psFree (view); 426 426 427 427 428 428 return true;
Note:
See TracChangeset
for help on using the changeset viewer.
