Changeset 21422 for trunk/psastro/src/psastroOneChipFit.c
- Timestamp:
- Feb 9, 2009, 11:25:34 AM (17 years ago)
- File:
-
- 1 edited
-
trunk/psastro/src/psastroOneChipFit.c (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psastro/src/psastroOneChipFit.c
r21409 r21422 6 6 * 7 7 * @author IfA 8 * @version $Revision: 1. 4$ $Name: not supported by cvs2svn $9 * @date $Date: 2009-02-0 7 02:03:34 $8 * @version $Revision: 1.5 $ $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 */ … … 23 23 24 24 // default value for match/fit : radius is in pixels 25 REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 25 REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 26 26 27 27 // run the match/fit sequence NITER times 28 REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 28 REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 29 29 30 30 // correct radius to FP units (physical pixel scale in microns per pixel) 31 REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 31 REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 32 32 RADIUS *= pixelScale; 33 33 … … 44 44 45 45 for (int iter = 0; iter < nIter; iter++) { 46 47 char name[128];48 46 49 sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter); 50 float radius = psMetadataLookupF32 (&status, recipe, name); 51 radius *= pixelScale; 52 if (!status || (radius == 0.0)) { 53 radius = RADIUS; 54 } 47 char name[128]; 48 49 sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter); 50 float radius = psMetadataLookupF32 (&status, recipe, name); 51 radius *= pixelScale; 52 if (!status || (radius == 0.0)) { 53 radius = RADIUS; 54 } 55 55 56 56 57 58 59 60 61 62 57 // use small radius to match stars 58 match = pmAstromRadiusMatchFP (rawstars, refstars, radius); 59 if (match == NULL) { 60 psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); 61 return false; 62 } 63 63 64 65 66 67 64 // modify the order to correspond to the actual number of matched stars: 65 int Ndof_min = 3; 66 int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1)); 67 order = PS_MIN (order, order_max); 68 68 69 70 71 69 // if ((match->n < 11) && (order >= 3)) order = 2; 70 // if ((match->n < 7) && (order >= 2)) order = 1; 71 // if ((match->n < 4) && (order >= 1)) order = 0; 72 72 73 74 psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 75 76 return false; 77 } 73 if (order < 1) { 74 psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 75 psFree (match); 76 return false; 77 } 78 78 79 80 81 82 83 84 85 86 87 88 89 79 // create output toFPA; set masks appropriate to the Elixir DVO astrometry format 80 psFree (chip->toFPA); 81 chip->toFPA = psPlaneTransformAlloc (order, order); 82 for (int i = 0; i <= chip->toFPA->x->nX; i++) { 83 for (int j = 0; j <= chip->toFPA->x->nY; j++) { 84 if (i + j > order) { 85 chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET; 86 chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; 87 } 88 } 89 } 90 90 91 92 93 94 95 91 // XXX allow statitic to be set by the user 92 // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); 93 fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); 94 fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA"); 95 fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER"); 96 96 97 // improved fit for astrometric terms 98 results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); 99 if (!results) { 100 psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); 101 psFree (match); 102 psFree (fitStats); 103 return false; 104 } 105 106 // determine fromFPA transformation and apply new transformation to raw & ref stars 107 psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); 108 109 // toSky converts from FPA & TPA units (microns) to sky units (radians) 110 float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; 111 // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; 112 float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; 113 int astNstar = results->yStats->clippedNvalues; 114 psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); 97 // improved fit for astrometric terms 98 results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); 99 if (!results) { 100 psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); 101 psFree (match); 102 psFree (fitStats); 103 return false; 104 } 115 105 116 if (iter < nIter - 1) { 117 psFree (fitStats); 118 psFree (results); 119 psFree (match); 120 } 106 // determine fromFPA transformation and apply new transformation to raw & ref stars 107 psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); 108 109 // toSky converts from FPA & TPA units (microns) to sky units (radians) 110 float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; 111 // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; 112 float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; 113 int astNstar = results->yStats->clippedNvalues; 114 psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); 115 116 if (iter < nIter - 1) { 117 psFree (fitStats); 118 psFree (results); 119 psFree (match); 120 } 121 121 } 122 122 … … 139 139 if (astError > maxError) { 140 140 psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); 141 141 validSolution = false; 142 142 } 143 143 if (astNstar < minNstar) { 144 144 psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); 145 145 validSolution = false; 146 146 } 147 147 … … 150 150 psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "astrometry error (arcsec)", astError); 151 151 if (validSolution) { 152 153 152 psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar)); 153 psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); 154 154 } else { 155 156 155 psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); 156 psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); 157 157 } 158 158 psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars … … 160 160 // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars 161 161 // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); 162 162 163 163 // XXX check if we correctly applied the new transformation: 164 164 if (psTraceGetLevel("psastro.dump") > 0) { 165 166 167 165 psastroDumpRawstars (rawstars, fpa, chip); 166 psastroDumpMatchedStars ("match.dat", rawstars, refstars, match); 167 psastroDumpStars (refstars, "refstars.cal.dat"); 168 168 } 169 169 170 p sastroVisualPlotOneChipFit (rawstars, refstars, match, recipe);170 pmAstromVisualPlotOneChipFit (rawstars, refstars, match, recipe); 171 171 172 172 if (psTraceGetLevel("psastro.plot") > 0) { 173 173 psastroPlotOneChipFit (rawstars, refstars, match, recipe); 174 174 } 175 175
Note:
See TracChangeset
for help on using the changeset viewer.
