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