- Timestamp:
- Oct 18, 2008, 2:26:45 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/cnb_branch_20081011/psastro/src/psastroAstromGuess.c
r19977 r20263 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 // dump or plot the resulting projected positions 141 if (psTraceGetLevel("psastro.dump") > 0) { 142 psastroDumpRawstars (rawstars, fpa, chip); 143 } 144 145 if (psTraceGetLevel("psastro.plot") > 0) { 146 psastroPlotRawstars (rawstars, fpa, chip, recipe); 147 } 140 // dump or plot the resulting projected positions 141 if (psTraceGetLevel("psastro.dump") > 0) { 142 psastroDumpRawstars (rawstars, fpa, chip); 143 } 144 145 psastroVisualPlotRawStars(rawstars, fpa, chip, recipe); 146 147 if (psTraceGetLevel("psastro.plot") > 0) { 148 psastroPlotRawstars (rawstars, fpa, chip, recipe); 149 } 148 150 } 149 151 } … … 155 157 psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR", PS_META_REPLACE, "", *nStars); 156 158 if (*nStars == 0) { 157 158 159 160 } 161 162 psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 163 DEG_RAD*RAmin, DEG_RAD*DECmin, 164 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); 165 167 166 168 psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", RAmin); … … 201 203 pmHDU *hdu = pmFPAviewThisHDU (view, fpa); 202 204 if (bilevelAstrometry) { 203 204 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 205 206 } 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 } 207 209 } else { 208 209 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 210 211 } 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 } 212 214 } 213 215 return true; … … 223 225 // load mosaic-level astrometry? 224 226 if (phu) { 225 226 227 228 227 char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); 228 if (ctype) { 229 *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); 230 } 229 231 } 230 232 if (*bilevelAstrometry) { 231 232 } 233 pmAstromReadBilevelMosaic (fpa, phu->header); 234 } 233 235 psFree (view); 234 236 return true; … … 243 245 pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); 244 246 if (!input) { 245 246 247 psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); 248 return false; 247 249 } 248 250 … … 273 275 if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } 274 276 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 277 psPlane ptCH, ptFP, ptTP; 278 psSphere ptSky; 279 280 ptCH.x = 0; 281 ptCH.y = 0; 282 psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); 283 psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); 284 psDeproject (&ptSky, &ptTP, fpa->toSky); 285 psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); 286 287 // new corner locations based on the calibrated astrometry 288 psVectorAppend (cornerLn, ptFP.x); 289 psVectorAppend (cornerMn, ptFP.y); 290 psVectorAppend (cornerPn, ptTP.x); 291 psVectorAppend (cornerQn, ptTP.y); 292 psVectorAppend (cornerRn, ptSky.r); 293 psVectorAppend (cornerDn, ptSky.d); 292 294 } 293 295 … … 298 300 299 301 for (int i = 0; i < cornerRo->n; i++) { 300 301 302 303 304 305 306 307 308 309 302 303 psPlane ptTP; 304 psSphere ptSky; 305 306 ptSky.r = cornerRo->data.F32[i]; 307 ptSky.d = cornerDo->data.F32[i]; 308 309 psProject (&ptTP, &ptSky, fpa->toSky); 310 psVectorAppend (cornerPs, ptTP.x); 311 psVectorAppend (cornerQs, ptTP.y); 310 312 } 311 313 … … 313 315 map->x->coeffMask[1][1] = PS_POLY_MASK_SET; 314 316 map->y->coeffMask[1][1] = PS_POLY_MASK_SET; 315 317 316 318 psVectorFitPolynomial2D (map->x, NULL, 0, cornerPn, NULL, cornerPs, cornerQs); 317 319 psVectorFitPolynomial2D (map->y, NULL, 0, cornerQn, NULL, cornerPs, cornerQs); 318 320 319 321 // apply the linear fit... 320 322 psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs); … … 325 327 psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf); 326 328 329 psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); 330 327 331 psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); 328 332 psStats *statsQ = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); … … 333 337 float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]); 334 338 float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]); 335 339 336 340 psLogMsg ("psastro", 3, "boresite offset : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]); 337 341 psLogMsg ("psastro", 3, "boresite angle : %f, scale: %f", angle*PS_DEG_RAD, scale); … … 341 345 psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER"); 342 346 if (!header) { 343 344 345 347 header = psMetadataAlloc(); 348 psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); 349 psFree (header); // drop this reference 346 350 } 347 351 … … 354 358 355 359 if (0) { 356 357 358 359 cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], 360 361 362 360 FILE *f = fopen ("corners.dat", "w"); 361 for (int i = 0; i < cornerRo->n; i++) { 362 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", 363 cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], 364 cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); 365 } 366 fclose (f); 363 367 } 364 368 … … 381 385 psFree (map); 382 386 psFree (view); 383 387 384 388 385 389 return true;
Note:
See TracChangeset
for help on using the changeset viewer.
