Changeset 21422 for trunk/psastro/src/psastroAstromGuess.c
- Timestamp:
- Feb 9, 2009, 11:25:34 AM (17 years ago)
- File:
-
- 1 edited
-
trunk/psastro/src/psastroAstromGuess.c (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psastro/src/psastroAstromGuess.c
r21409 r21422 6 6 * 7 7 * @author IfA 8 * @version $Revision: 1.3 4$ $Name: not supported by cvs2svn $9 * @date $Date: 2009-02-0 7 02:03:34 $8 * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2009-02-09 21:25:34 $ 10 10 * Copyright 2009 Institute for Astronomy, University of Hawaii 11 11 */ … … 42 42 psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE); 43 43 if (!recipe) { 44 45 44 psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); 45 return false; 46 46 } 47 47 … … 49 49 bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); 50 50 if (!status) { 51 51 useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); 52 52 } 53 53 … … 55 55 pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); 56 56 if (!input) { 57 58 57 psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); 58 return false; 59 59 } 60 60 … … 62 62 double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); 63 63 if (!status) { 64 psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); 65 return false; 66 } 64 psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); 65 return false; 66 } 67 67 68 68 psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32); … … 80 80 bool bilevelAstrometry = false; 81 81 if (!useModel) { 82 82 psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); 83 83 } 84 84 … … 88 88 if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } 89 89 90 91 92 90 if (!useModel) { 91 if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; 92 } 93 93 94 94 if (newFPA) { 95 95 newFPA = false; 96 97 96 while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; 97 while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; 98 98 RAminSky = fpa->toSky->R - M_PI; 99 99 RAmaxSky = fpa->toSky->R + M_PI; 100 100 } 101 101 102 103 { 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 102 // report and save the current best guess for the chip 0,0 pixel coordinates 103 { 104 psPlane ptCH, ptFP, ptTP; 105 psSphere ptSky; 106 107 ptCH.x = 0; 108 ptCH.y = 0; 109 psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); 110 psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); 111 psDeproject (&ptSky, &ptTP, fpa->toSky); 112 psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); 113 114 psVectorAppend (cornerL, ptFP.x); 115 psVectorAppend (cornerM, ptFP.y); 116 psVectorAppend (cornerP, ptTP.x); 117 psVectorAppend (cornerQ, ptTP.y); 118 psVectorAppend (cornerR, ptSky.r); 119 psVectorAppend (cornerD, ptSky.d); 120 } 121 121 122 122 // apply the new WCS guess data to all of the data in the readouts … … 132 132 if (rawstars == NULL) { continue; } 133 133 134 134 *nStars += rawstars->n; 135 135 for (int i = 0; i < rawstars->n; i++) { 136 136 pmAstromObj *raw = rawstars->data[i]; … … 151 151 } 152 152 153 154 155 156 157 158 p sastroVisualPlotRawStars(rawstars, fpa, chip, recipe);159 160 161 162 153 // dump or plot the resulting projected positions 154 if (psTraceGetLevel("psastro.dump") > 0) { 155 psastroDumpRawstars (rawstars, fpa, chip); 156 } 157 158 pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe); 159 160 if (psTraceGetLevel("psastro.plot") > 0) { 161 psastroPlotRawstars (rawstars, fpa, chip, recipe); 162 } 163 163 } 164 164 } … … 170 170 psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR", PS_META_REPLACE, "", *nStars); 171 171 if (*nStars == 0) { 172 173 174 175 } 176 177 psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 178 DEG_RAD*RAmin, DEG_RAD*DECmin, 179 172 psLogMsg ("psastro", 2, "no sources available for astrometry\n"); 173 psFree (view); 174 return true; 175 } 176 177 psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 178 DEG_RAD*RAmin, DEG_RAD*DECmin, 179 DEG_RAD*RAmax, DEG_RAD*DECmax); 180 180 181 181 psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", RAmin); … … 216 216 pmHDU *hdu = pmFPAviewThisHDU (view, fpa); 217 217 if (bilevelAstrometry) { 218 219 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 220 221 } 218 if (!pmAstromReadBilevelChip (chip, hdu->header)) { 219 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 220 return false; 221 } 222 222 } else { 223 224 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 225 226 } 223 if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { 224 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 225 return false; 226 } 227 227 } 228 228 return true; … … 238 238 // load mosaic-level astrometry? 239 239 if (phu) { 240 241 242 243 240 char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); 241 if (ctype) { 242 *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); 243 } 244 244 } 245 245 if (*bilevelAstrometry) { 246 247 } 246 pmAstromReadBilevelMosaic (fpa, phu->header); 247 } 248 248 psFree (view); 249 249 return true; … … 258 258 pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); 259 259 if (!input) { 260 261 260 psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); 261 return false; 262 262 } 263 263 … … 290 290 if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } 291 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 320 321 322 323 324 325 326 327 328 329 330 331 332 292 // XXX we are currently inconsistent with marking the good vs the bad data 293 // psastroChipAstrom sets data_exists to false if the fit is bad. this is 294 // probably wrong since it implies there is no data! 295 296 // skip chips for which the astrometry failed (NASTRO == 0) 297 if (!chip->cells->n) goto skip_chip; 298 pmCell *cell = chip->cells->data[0]; 299 if (!cell) goto skip_chip; 300 301 if (!cell->readouts->n) goto skip_chip; 302 pmReadout *readout = cell->readouts->data[0]; 303 if (!readout) goto skip_chip; 304 305 psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); 306 if (!updates) goto skip_chip; 307 308 int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); 309 if (!nAstro) goto skip_chip; 310 311 float astError = psMetadataLookupF32 (&status, updates, "CERROR"); 312 if (fabs(astError) < 1e-6) goto skip_chip; 313 314 psPlane ptCH, ptFP, ptTP; 315 psSphere ptSky; 316 317 ptCH.x = 0; 318 ptCH.y = 0; 319 psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); 320 psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); 321 psDeproject (&ptSky, &ptTP, fpa->toSky); 322 psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); 323 324 // new corner locations based on the calibrated astrometry 325 psVectorAppend (cornerLn, ptFP.x); 326 psVectorAppend (cornerMn, ptFP.y); 327 psVectorAppend (cornerPn, ptTP.x); 328 psVectorAppend (cornerQn, ptTP.y); 329 psVectorAppend (cornerRn, ptSky.r); 330 psVectorAppend (cornerDn, ptSky.d); 331 psVectorAppend (cornerMK, 0); 332 continue; 333 333 334 334 skip_chip: 335 336 337 338 339 340 341 342 335 // new corner locations based on the calibrated astrometry 336 psVectorAppend (cornerLn, 0.0); 337 psVectorAppend (cornerMn, 0.0); 338 psVectorAppend (cornerPn, 0.0); 339 psVectorAppend (cornerQn, 0.0); 340 psVectorAppend (cornerRn, 0.0); 341 psVectorAppend (cornerDn, 0.0); 342 psVectorAppend (cornerMK, 1); 343 343 } 344 344 … … 349 349 350 350 for (int i = 0; i < cornerRo->n; i++) { 351 352 353 354 355 356 357 358 359 360 351 352 psPlane ptTP; 353 psSphere ptSky; 354 355 ptSky.r = cornerRo->data.F32[i]; 356 ptSky.d = cornerDo->data.F32[i]; 357 358 psProject (&ptTP, &ptSky, fpa->toSky); 359 psVectorAppend (cornerPs, ptTP.x); 360 psVectorAppend (cornerQs, ptTP.y); 361 361 } 362 362 … … 364 364 map->x->coeffMask[1][1] = PS_POLY_MASK_SET; 365 365 map->y->coeffMask[1][1] = PS_POLY_MASK_SET; 366 366 367 367 // fit the valid chips, mask the invalid chips 368 368 psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs); 369 369 psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs); 370 370 371 371 // apply the linear fit... 372 372 psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs); … … 377 377 psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf); 378 378 379 p sastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);379 pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); 380 380 381 381 psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); … … 387 387 float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]); 388 388 float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]); 389 389 390 390 psLogMsg ("psastro", 3, "boresite offset : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]); 391 391 psLogMsg ("psastro", 3, "boresite angle : %f, scale: %f", angle*PS_DEG_RAD, scale); … … 395 395 psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER"); 396 396 if (!header) { 397 398 399 397 header = psMetadataAlloc(); 398 psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); 399 psFree (header); // drop this reference 400 400 } 401 401 … … 408 408 409 409 if (DEBUG) { 410 411 412 413 cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], 414 415 416 410 FILE *f = fopen ("corners.dat", "w"); 411 for (int i = 0; i < cornerRo->n; i++) { 412 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", 413 cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], 414 cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); 415 } 416 fclose (f); 417 417 } 418 418 … … 437 437 psFree (map); 438 438 psFree (view); 439 439 440 440 441 441 return true;
Note:
See TracChangeset
for help on using the changeset viewer.
