Changeset 21422 for trunk/psastro/src/psastroFixChips.c
- Timestamp:
- Feb 9, 2009, 11:25:34 AM (17 years ago)
- File:
-
- 1 edited
-
trunk/psastro/src/psastroFixChips.c (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psastro/src/psastroFixChips.c
r21409 r21422 1 1 /** @file psastroFixChips.c 2 2 * 3 * @brief 3 * @brief 4 4 * 5 5 * @ingroup libpsastro 6 6 * 7 7 * @author IfA 8 * @version $Revision: 1. 8$ $Name: not supported by cvs2svn $9 * @date $Date: 2009-02-0 7 02:03:34 $8 * @version $Revision: 1.9 $ $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 */ … … 25 25 bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS"); 26 26 if (!status) { 27 27 fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); 28 28 } 29 29 if (!fixChips) return true; … … 37 37 pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); 38 38 if (!input) { 39 40 39 psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); 40 return false; 41 41 } 42 42 … … 60 60 // files associated with the science image 61 61 if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { 62 63 62 psError (PS_ERR_IO, false, "Can't load the astrometry model file"); 63 return false; 64 64 } 65 65 … … 76 76 77 77 if (DEBUG) { 78 79 78 f = fopen ("corners.raw.dat", "w"); 79 chipName = NULL; 80 80 } 81 81 82 82 pmChip *obsChip = NULL; 83 83 while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 84 if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } 85 86 // XXX we are currently inconsistent with marking the good vs the bad data 87 // psastroChipAstrom sets data_exists to false if the fit is bad. this is 88 // probably wrong since it implies there is no data! 89 90 // skip chips for which the astrometry failed (NASTRO == 0) 91 if (!obsChip->cells->n) continue; 92 pmCell *cell = obsChip->cells->data[0]; 93 if (!cell) continue; 94 95 if (!cell->readouts->n) continue; 96 pmReadout *readout = cell->readouts->data[0]; 97 if (!readout) continue; 98 99 psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); 100 if (!updates) continue; 101 102 int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); 103 if (!nAstro) continue; 104 105 // set the chip astrometry using the astrom file 106 pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); 107 108 psRegion *region = pmChipPixels (obsChip); 109 psPlane ptCP, ptFP; 110 111 ptCP.x = region->x0; ptCP.y = region->y0; 112 psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); 113 xObs->data.F32[nPts] = ptFP.x; 114 yObs->data.F32[nPts] = ptFP.y; 115 psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); 116 xRef->data.F32[nPts] = ptFP.x; 117 yRef->data.F32[nPts] = ptFP.y; 118 119 if (DEBUG) { 120 chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); 121 fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); 122 } 123 nPts ++; 124 125 ptCP.x = region->x0; ptCP.y = region->y1; 126 psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); 127 xObs->data.F32[nPts] = ptFP.x; 128 yObs->data.F32[nPts] = ptFP.y; 129 psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); 130 xRef->data.F32[nPts] = ptFP.x; 131 yRef->data.F32[nPts] = ptFP.y; 132 133 if (DEBUG) { 134 chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); 135 fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); 136 } 137 nPts ++; 138 139 ptCP.x = region->x1; ptCP.y = region->y1; 140 psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); 141 xObs->data.F32[nPts] = ptFP.x; 142 yObs->data.F32[nPts] = ptFP.y; 143 psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); 144 xRef->data.F32[nPts] = ptFP.x; 145 yRef->data.F32[nPts] = ptFP.y; 146 147 if (DEBUG) { 148 chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); 149 fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); 150 } 151 nPts ++; 152 153 ptCP.x = region->x1; ptCP.y = region->y0; 154 psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP); 155 xObs->data.F32[nPts] = ptFP.x; 156 yObs->data.F32[nPts] = ptFP.y; 157 psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP); 158 xRef->data.F32[nPts] = ptFP.x; 159 yRef->data.F32[nPts] = ptFP.y; 160 161 if (DEBUG) { 162 chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME"); 163 fprintf (f, "%s %f %f %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]); 164 } 165 nPts ++; 166 167 psFree (region); 168 168 } 169 169 xObs->n = yObs->n = xRef->n = yRef->n = nPts; 170 170 if (DEBUG) fclose (f); 171 171 172 172 psPlaneTransform *map = psPlaneTransformAlloc (1, 1); 173 173 174 174 psVector *mask = psVectorAlloc (nPts, PS_TYPE_VECTOR_MASK); 175 175 psVectorInit (mask, 0); … … 179 179 180 180 for (int i = 0; i < 3; i++) { 181 182 183 184 185 181 psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef); 182 psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n); 183 184 psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef); 185 psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n); 186 186 } 187 187 188 188 // loop over all chips, select the outliers, and replace the measured astrometry with the model 189 // the measured transformation above must be applied to make the comparison, and also then applied to the 189 // the measured transformation above must be applied to make the comparison, and also then applied to the 190 190 // model transformation 191 191 192 192 if (DEBUG) { 193 194 195 196 197 198 199 200 201 193 f = fopen ("corners.fit.dat", "w"); 194 for (int i = 0; i < xObs->n; i++) { 195 psPlane obsCoord, refCoord; 196 refCoord.x = xRef->data.F32[i]; 197 refCoord.y = yRef->data.F32[i]; 198 psPlaneTransformApply (&obsCoord, map, &refCoord); 199 fprintf (f, "%f %f %f %f %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y); 200 } 201 fclose (f); 202 202 } 203 203 … … 211 211 212 212 while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 // XXX for now, let's just fail on the bad chips. In the future, let's try to recover, but we still need to 272 273 274 275 276 277 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 psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", 317 318 319 } 320 321 322 323 324 325 326 327 328 329 330 331 } 332 333 p sastroVisualPlotFixChips (input, xObs, yObs);213 psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); 214 if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } 215 216 // set the chip astrometry using the astrom file 217 pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); 218 219 // bad Astrometry test: ref pixel or angle outside nominal 220 221 psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; 222 psPlane obsCoord, refCoord, tmpCoord; 223 224 // find location of 0,0 pixel in focal plane coords for this chip 225 psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); 226 227 // find location of 0,0 pixel in focal plane coords for ref chip 228 // apply the global field rotation and offset before comparing 229 psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel); 230 psPlaneTransformApply (&refCoord, map, &tmpCoord); 231 232 psPlane offPixel = {100.0, 0.0, 100.0, 0.0}; 233 psPlane obsOffPt, refOffPt; 234 235 // find location of 0,0 pixel in focal plane coords for this chip 236 psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); 237 238 // find location of 0,0 pixel in focal plane coords for ref chip 239 psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel); 240 psPlaneTransformApply (&refOffPt, map, &tmpCoord); 241 242 double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); 243 double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x); 244 245 bool badAstrom = false; 246 badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; 247 badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; 248 badAstrom |= fabs(obsAngle - refAngle) > angleTol; 249 250 fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); 251 252 // XXX for now, just use first readout 253 pmCell *cell = obsChip->cells->data[0]; 254 pmReadout *readout = cell->readouts->data[0]; 255 256 // update the header (pull or create local view to entry on readout->analysis) 257 psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); 258 if (!updates) { 259 updates = psMetadataAlloc (); 260 psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); 261 psFree (updates); 262 } 263 264 psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x); 265 psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y); 266 psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle); 267 268 // for successful chips, save the measured offsets in the header 269 if (!badAstrom) continue; 270 271 // XXX for now, let's just fail on the bad chips. In the future, let's try to recover, but we still need to 272 // catch the failures relative to the model 273 psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); 274 continue; 275 276 psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n", 277 view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); 278 279 psFree (obsChip->toFPA); 280 psFree (obsChip->fromFPA); 281 282 // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter 283 // XXX this only works if toTPA is at most a linear transformation 284 psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY); 285 for (int i = 0; i <= refChip->toFPA->x->nX; i++) { 286 for (int j = 0; j <= refChip->toFPA->x->nY; j++) { 287 double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; 288 double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; 289 toFPA->x->coeff[i][j] = f1 + f2; 290 291 double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j]; 292 double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j]; 293 toFPA->y->coeff[i][j] = g1 + g2; 294 } 295 } 296 toFPA->x->coeff[0][0] += map->x->coeff[0][0]; 297 toFPA->y->coeff[0][0] += map->y->coeff[0][0]; 298 299 psRegion *region = pmChipPixels (obsChip); 300 obsChip->toFPA = toFPA; 301 obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50); 302 psFree (region); 303 304 // use the new position to re-try the match fit 305 // select the raw objects for this readout 306 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); 307 if (rawstars == NULL) { continue; } 308 309 // select the raw objects for this readout 310 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); 311 if (refstars == NULL) { continue; } 312 313 // the absolute minimum number of stars is 4 (for order = 1) 314 if ((rawstars->n < 4) || (refstars->n < 4)) { 315 readout->data_exists = false; 316 psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", 317 rawstars->n, refstars->n); 318 continue; 319 } 320 321 psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars); 322 323 // XXX update the header with info to reflect the failure 324 if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) { 325 readout->data_exists = false; 326 psLogMsg ("psastro", 3, "failed to find a solution\n"); 327 continue; 328 } 329 330 pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL); 331 } 332 333 pmAstromVisualPlotFixChips (input, xObs, yObs); 334 334 psFree (xObs); 335 335 psFree (yObs);
Note:
See TracChangeset
for help on using the changeset viewer.
