Index: /trunk/psastro/src/Makefile.am =================================================================== --- /trunk/psastro/src/Makefile.am (revision 10829) +++ /trunk/psastro/src/Makefile.am (revision 10830) @@ -12,5 +12,7 @@ libpsastro_la_SOURCES = \ -psastroArguments.c \ +psastroArguments.c \ +psastroErrorCodes.c \ +psastroVersion.c \ psastroCleanup.c \ psastroParseCamera.c \ @@ -18,6 +20,6 @@ psastroDataSave.c \ psastroAstromGuess.c \ -psastroLoadRefstars.c \ -psastroChooseRefstars.c \ +psastroLoadRefstars.c \ +psastroChooseRefstars.c \ psastroConvert.c \ psastroChipAstrom.c \ @@ -28,12 +30,9 @@ psastroRefstarSubset.c \ psastroMosaicAstrom.c \ -psastroMosaicGetGrads.c \ +psastroMosaicGradients.c \ psastroMosaicChipAstrom.c \ +psastroMosaicOneChip.c \ psastroMosaicSetAstrom.c \ -psastroMosaicSetMatch.c \ -psastroMosaicHeaders.c \ -psastroMosaicRescaleChips.c \ -psastroErrorCodes.c \ -psastroVersion.c +psastroMosaicSetMatch.c include_HEADERS = \ Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 10829) +++ /trunk/psastro/src/psastro.h (revision 10830) @@ -56,16 +56,16 @@ // mosaic fitting functions -psArray *psastroMosaicGetGrads (pmFPA *fpa, psMetadata *recipe); +psArray *psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe); +bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe); bool psastroMosaicAstrom (pmConfig *config, psArray *refs); -bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe); -bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe); +bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration); +bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration); bool psastroMosaicSetAstrom (pmFPA *fpa); bool psastroMosaicHeaders (pmConfig *config); bool psastroMosaicRescaleChips (pmFPA *fpa); - -bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale); -psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale); +bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration); // Return version strings. psString psastroVersion(void); psString psastroVersionLong(void); + Index: /trunk/psastro/src/psastroAstromGuess.c =================================================================== --- /trunk/psastro/src/psastroAstromGuess.c (revision 10829) +++ /trunk/psastro/src/psastroAstromGuess.c (revision 10830) @@ -11,5 +11,4 @@ bool newFPA = true; bool status = false; - bool isMosaic = false; double RAmin = NAN; double RAmax = NAN; @@ -25,7 +24,4 @@ DECmin = DECmax = -90; RAmin = RAmax = RAminSky = RAmaxSky = 0; - - // is this a mosaic astrometry analysis? - psMetadataLookupStr (&isMosaic, config->arguments, "MOSASTRO"); // select the current recipe @@ -43,7 +39,10 @@ } - double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE"); - if (!status) plateScale = 1.0; - plateScale = 1.0; + // physical pixel scale in microns per pixel + double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); + if (!status) { + psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); + return false; + } pmFPAview *view = pmFPAviewAlloc (0); @@ -57,5 +56,5 @@ pmHDU *hdu = pmFPAviewThisHDU (view, fpa); - pmAstromReadWCS (fpa, chip, hdu->header, plateScale, isMosaic); + pmAstromReadWCS (fpa, chip, hdu->header, pixelScale); if (newFPA) { newFPA = false; @@ -69,5 +68,4 @@ // apply the new WCS guess data to all of the data in the readouts - // XXX should this go into a different function? this would separate WCS interpretation from application while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); @@ -85,5 +83,5 @@ psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); - psPlaneDistortApply (raw->TP, fpa->toTPA, raw->FP, 0.0, 0.0); + psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP); psDeproject (raw->sky, raw->TP, fpa->toSky); @@ -100,5 +98,5 @@ psProject (tp, raw->sky, fpa->toSky); - psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0); + psPlaneTransformApply (fp, fpa->fromTPA, tp); psPlaneTransformApply (ch, chip->fromFPA, fp); Index: /trunk/psastro/src/psastroChooseRefstars.c =================================================================== --- /trunk/psastro/src/psastroChooseRefstars.c (revision 10829) +++ /trunk/psastro/src/psastroChooseRefstars.c (revision 10830) @@ -62,5 +62,5 @@ psProject (ref->TP, ref->sky, fpa->toSky); - psPlaneDistortApply (ref->FP, fpa->fromTPA, ref->TP, 0.0, 0.0); + psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP); psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); Index: /trunk/psastro/src/psastroMosaicAstrom.c =================================================================== --- /trunk/psastro/src/psastroMosaicAstrom.c (revision 10829) +++ /trunk/psastro/src/psastroMosaicAstrom.c (revision 10830) @@ -1,8 +1,10 @@ # include "psastro.h" +# define NONLIN_TOL 0.001 /* tolerance in pixels */ +// XXX require this fpa to have multiple chip extensions and a PHU? bool psastroMosaicAstrom (pmConfig *config, psArray *refs) { bool status; - psArray *grads = NULL; + psArray *gradients = NULL; // select the current recipe @@ -20,44 +22,71 @@ } - int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.DISTORTION.ORDER"); - if (!status) order = 3.0; - pmFPA *fpa = input->fpa; - psastroMosaicSetMatch (fpa, recipe); + // XXX before we do object matches, we need to fix failed chips + // compare chips with supplied mosaic model + // adjust significant outliers to match model - // XXX no input distortion model yet; existing fpa distortion is identity: replace it - // XXX make this a test if these are NULL; the load step can NULL them if nothing is loaded - // XXX if we have an input distortion model, the gradient is measuring the error rel to that model - // XXX how does this couple to the individual chip fits? + // given the existing per-chip astrometry, determine matches between raw and ref stars + psastroMosaicSetMatch (fpa, recipe, 0); + + // fitted chips will follow the local plate-scale, hiding the distortion + // modify the chip->toFPA scaling to match knowledge about pixel scale + psastroMosaicCommonScale (fpa, recipe); + + gradients = psastroMosaicGradients (fpa, recipe); + + // allocate mosaic-level polynomial transformation and set masks needed by DVO + int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER"); + if (!status) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n"); + return false; + } psFree (fpa->toTPA); - psFree (fpa->fromTPA); - fpa->toTPA = psPlaneDistortIdentity (order); - fpa->fromTPA = psPlaneDistortIdentity (order); + fpa->toTPA = psPlaneTransformAlloc (order, order); + for (int i = 0; i <= fpa->toTPA->x->nX; i++) { + for (int j = 0; j <= fpa->toTPA->x->nY; j++) { + if (i + j > order) { + fpa->toTPA->x->mask[i][j] = 1; + fpa->toTPA->y->mask[i][j] = 1; + } + } + } - grads = psastroMosaicGetGrads (fpa, recipe); - - // fit the measured gradients with the telescope distortion model (3rd order polynomial) - pmAstromFitDistortion (fpa, grads, recipe); - psastroMosaicRescaleChips (fpa); + // fit the measured gradients with the telescope distortion model (polynomial order based on toTPA) + pmAstromFitDistortion (fpa, gradients, recipe); + psFree (gradients); // apply the new distortion terms up and down // re-perform the match with a slightly tighter circle + // XXX modify match radius + // XXX set chip.order to 1 psastroMosaicSetAstrom (fpa); - psastroMosaicSetMatch (fpa, recipe); + psastroMosaicSetMatch (fpa, recipe, 1); + psastroMosaicChipAstrom (fpa, recipe, 1); - // fit the chips with linear fits - // re-match - psastroMosaicChipAstrom (fpa, recipe); - psastroMosaicSetMatch (fpa, recipe); + // do a second pass on the distortion with improved chip positions and rotations + // XXX do we need to iterate with this step until the distortions don't change? + psastroMosaicCommonScale (fpa, recipe); + gradients = psastroMosaicGradients (fpa, recipe); + pmAstromFitDistortion (fpa, gradients, recipe); + psFree (gradients); + + // now fit the chips under the common distortion with higher-order terms + // XXX modify match radius + // XXX set chip.order to NORDER + psastroMosaicSetAstrom (fpa); + psastroMosaicSetMatch (fpa, recipe, 2); + psastroMosaicChipAstrom (fpa, recipe, 2); - // re-fit the chips with higher-order fits - // re-match - psastroMosaicChipAstrom (fpa, recipe); - psastroMosaicSetMatch (fpa, recipe); + // save WCS and analysis metadata in update header + // XXX need to add global summary statistics + psMetadata *updates = psMetadataAlloc(); + pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL); + psMetadataAdd (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_DATA_METADATA, "psastro header stats", updates); + psFree (updates); - // re-fit the chips with higher-order fits - psastroMosaicChipAstrom (fpa, recipe); - psastroMosaicHeaders (config); + // update the headers based on the results + // psastroMosaicHeaders (config); return true; @@ -65,9 +94,9 @@ /* coordinate frame hierachy - pixels (on a given readout) - cell - chip - FP (focal plane) - TP (tangent plane) - sky (ra, dec) -*/ + * pixels (on a given readout) + * cell + * chip + * FP (focal plane) + * TP (tangent plane) + * sky (ra, dec) + */ Index: /trunk/psastro/src/psastroMosaicChipAstrom.c =================================================================== --- /trunk/psastro/src/psastroMosaicChipAstrom.c (revision 10829) +++ /trunk/psastro/src/psastroMosaicChipAstrom.c (revision 10830) @@ -1,5 +1,6 @@ # include "psastro.h" +# define NONLIN_TOL 0.001 /* tolerance in pixels */ -bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe) { +bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration) { pmChip *chip = NULL; @@ -22,16 +23,12 @@ if (! readout->data_exists) { continue; } - // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); - if (rawstars == NULL) { continue; } + // save WCS and analysis metadata in update header + psMetadata *updates = psMetadataAlloc(); - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); - if (refstars == NULL) { continue; } + psastroMosaicOneChip (chip, readout, recipe, updates, iteration); - psArray *match = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); - if (match == NULL) { continue; } - - // need to pass in an update header, sent in from above - pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, recipe, NULL); + pmAstromWriteBilevelChip (updates, chip, NONLIN_TOL); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_DATA_METADATA, "psastro header stats", updates); + psFree (updates); } } @@ -39,2 +36,5 @@ return true; } + +/* the iteration value may be 1 or 2. + */ Index: /trunk/psastro/src/psastroMosaicGradients.c =================================================================== --- /trunk/psastro/src/psastroMosaicGradients.c (revision 10830) +++ /trunk/psastro/src/psastroMosaicGradients.c (revision 10830) @@ -0,0 +1,44 @@ +# include "psastro.h" + +psArray *psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe) { + + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + psArray *gradients = NULL; + + pmFPAview *view = pmFPAviewAlloc (0); + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + // XXX there can only be one readout per chip, right? + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) { continue; } + + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + + psArray *match = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); + if (match == NULL) { continue; } + + // measure the local gradients for this set of stars + // XXX update the function prototype to accept an incoming gradient structure to which the new elements are added + gradients = pmAstromMeasureGradients (gradients, rawstars, refstars, match, recipe); + } + } + } + return (gradients); +} Index: /trunk/psastro/src/psastroMosaicOneChip.c =================================================================== --- /trunk/psastro/src/psastroMosaicOneChip.c (revision 10830) +++ /trunk/psastro/src/psastroMosaicOneChip.c (revision 10830) @@ -0,0 +1,103 @@ +# include "psastro.h" + +bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration) { + + bool status; + + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); + PS_ASSERT_PTR_NON_NULL(updates, false); + + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) return false; + + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) return false; + + psArray *match = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); + if (match == NULL) return false; + + if ((iteration < 1) || (iteration > 2)) { + psError(PSASTRO_ERR_UNKNOWN, false, "invalid iteration number %d\n", iteration); + return false; + } + + // correct to FP units (physical pixel scale in microns per pixel) + double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); + if (!status) { + psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); + return false; + } + + // set the order of the per-chip fit (iteration 1: 1, iteration 2: from recipe) + int order = 1; + if (iteration == 2) { + order = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.ORDER"); + if (!status) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic chip-level fit order\n"); + return false; + } + } + psFree (chip->toFPA); + chip->toFPA = psPlaneTransformAlloc (order, order); + for (int i = 0; i <= chip->toFPA->x->nX; i++) { + for (int j = 0; j <= chip->toFPA->x->nY; j++) { + if (i + j > order) { + chip->toFPA->x->mask[i][j] = 1; + chip->toFPA->y->mask[i][j] = 1; + } + } + } + + // XXX allow statitic to be set by the user + psStats *fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); + fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); + + // need to pass in an update header, sent in from above + pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); + if (!results) { + psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n"); + return false; + } + + // allowed limits for valid solutions + // XXX have these values depend on the iteration? + float maxError = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.MAX.ERROR"); + int minNstar = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.MIN.NSTAR"); + + // astError is the average 1D scatter in pixels ('results' are in FPA units = microns) + // XXX currently using supplied pixel scale: use toFPA->x->coeff[1][0], etc? + float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale; + int astNstar = results->yStats->clippedNvalues; + + bool validSolution = true; + + // XXX should these result in errors or be handled another way? + psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f, Nstars: %d", astError, astNstar); + if (astError > maxError) { + psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); + validSolution = false; + } + if (astNstar < minNstar) { + psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); + validSolution = false; + } + + // DVO expects NASTRO = 0 if we fail to find a solution + if (validSolution) { + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", astNstar); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(astNstar)); + } else { + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", 0); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", 0.0); + } + psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "", astError); + psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars + + psFree (fitStats); + psFree (results); + return true; +} Index: /trunk/psastro/src/psastroMosaicSetAstrom.c =================================================================== --- /trunk/psastro/src/psastroMosaicSetAstrom.c (revision 10829) +++ /trunk/psastro/src/psastroMosaicSetAstrom.c (revision 10830) @@ -30,5 +30,5 @@ psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); - psPlaneDistortApply (raw->TP, fpa->toTPA, raw->FP, 0.0, 0.0); + psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP); psDeproject (raw->sky, raw->TP, fpa->toSky); } @@ -41,5 +41,5 @@ psProject (ref->TP, ref->sky, fpa->toSky); - psPlaneDistortApply (ref->FP, fpa->fromTPA, ref->TP, 0.0, 0.0); + psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP); psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); } Index: /trunk/psastro/src/psastroMosaicSetMatch.c =================================================================== --- /trunk/psastro/src/psastroMosaicSetMatch.c (revision 10829) +++ /trunk/psastro/src/psastroMosaicSetMatch.c (revision 10830) @@ -1,5 +1,5 @@ # include "psastro.h" -bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe) { +bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration) { pmChip *chip = NULL; @@ -7,10 +7,17 @@ pmReadout *readout = NULL; pmFPAview *view = pmFPAviewAlloc (0); - - FILE *f; - char name[128]; + char radiusWord[64]; FILE *g1 = fopen ("raw.ps.dat", "w"); FILE *g2 = fopen ("ref.ps.dat", "w"); + + // use small radius to match stars (assume starting astrometry is good) + bool status = false; + sprintf (radiusWord, "PSASTRO.MOSAIC.RADIUS.N%d", iteration); + double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord); + if (!status) { + psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord); + return NULL; + } // this loop selects the matched stars for all chips @@ -37,24 +44,6 @@ psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); - // use small radius to match stars (assume starting astrometry is good) - // XXX should this take a (double radius)? - psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, recipe); + psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS); psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n); - - sprintf (name, "raw.%02d.dat", view->chip); - f = fopen (name, "w"); - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - fprintf (f, "%f %f %f %f %f %f\n", raw->chip->x, raw->chip->y, raw->FP->x, raw->FP->y, raw->sky->r, raw->sky->d); - } - fclose (f); - - sprintf (name, "ref.%02d.dat", view->chip); - f = fopen (name, "w"); - for (int i = 0; i < refstars->n; i++) { - pmAstromObj *ref = refstars->data[i]; - fprintf (f, "%f %f %f %f %f %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->sky->r, ref->sky->d); - } - fclose (f); for (int i = 0; i < matches->n; i++) { Index: /trunk/psastro/src/psastroOneChip.c =================================================================== --- /trunk/psastro/src/psastroOneChip.c (revision 10829) +++ /trunk/psastro/src/psastroOneChip.c (revision 10830) @@ -25,6 +25,21 @@ psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + // supplied radius is in pixels + double RADIUS = psMetadataLookupF32 (&status, recipe, "PSASTRO.MATCH.RADIUS"); + if (!status) { + psError(PS_ERR_IO, false, "Failed to lookup matching radius"); + return false; + } + + // correct to FP units (physical pixel scale in microns per pixel) + double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); + if (!status) { + psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); + return false; + } + RADIUS *= pixelScale; + // use small radius to match stars - psArray *match = pmAstromRadiusMatch (rawstars, refstars, recipe); + psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS); if (stats == NULL) { psError(PSASTRO_ERR_UNKNOWN, false, "failed to find radius-matched sources\n"); @@ -32,25 +47,72 @@ } - int nX = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NX"); - if (!status) nX = 1; - - int nY = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NY"); - if (!status) nY = 1; - - psPlaneTransform *toFPAout = psPlaneTransformAlloc (nX, nY); - - // improved fit for astrometric terms - if (!pmAstromMatchFit (toFPAout, rawstars, refstars, match, recipe, updates)) { - psError(PSASTRO_ERR_DATA, false, "failed to find a valid fitted match solution\n"); + // create the output fit model + int order = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.ORDER"); + if (!status) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n"); return false; } psFree (chip->toFPA); - chip->toFPA = toFPAout; + chip->toFPA = psPlaneTransformAlloc (order, order); + for (int i = 0; i <= chip->toFPA->x->nX; i++) { + for (int j = 0; j <= chip->toFPA->x->nY; j++) { + if (i + j > order) { + chip->toFPA->x->mask[i][j] = 1; + chip->toFPA->y->mask[i][j] = 1; + } + } + } - // write results + // XXX allow statitic to be set by the user + psStats *fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA"); + fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER"); + + // improved fit for astrometric terms + pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); + if (!results) { + psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n"); + return false; + } + + // allowed limits for valid solutions + float maxError = psMetadataLookupF32 (&status, recipe, "PSASTRO.MAX.ERROR"); + int minNstar = psMetadataLookupS32 (&status, recipe, "PSASTRO.MIN.NSTAR"); + + // astError is the average 1D scatter in pixels ('results' are in FPA units = microns) + float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale; + int astNstar = results->yStats->clippedNvalues; + + bool validSolution = true; + + // XXX should these result in errors or be handled another way? + psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f, Nstars: %d", astError, astNstar); + if (astError > maxError) { + psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError); + validSolution = false; + } + if (astNstar < minNstar) { + psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar); + validSolution = false; + } + + // DVO expects NASTRO = 0 if we fail to find a solution + if (validSolution) { + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", astNstar); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(astNstar)); + } else { + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", 0); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", 0.0); + } + psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "", astError); + psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars + +// write results psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); psFree (match); psFree (stats); + psFree (results); + psFree (fitStats); psFree (gridStats); return true; Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 10829) +++ /trunk/psastro/src/psastroUtils.c (revision 10830) @@ -1,4 +1,76 @@ # include "psastro.h" # define RENORM 0 + +// fix this to look up the value in the chip concepts +static double getChipPixelScale (pmChip *chip) { + return 10.0; +} + +// I have an FPA structure with multiple chips. we have loaded or measured astrometry for each +// chip with independent pixel/degree scaling values. These will naturally compensate locally +// somewhat for the telescope distortion. to measure the telescope distortion, we need to +// force the chips to have the same pixel scale and measure the difference from that solution. +// Convert an FPA with disparate pixel scales to a common pixel scale (perhaps depending on the +// chip -- eg, TC3) + +bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe) { + + // options : use the MIN or MAX chip as global reference or supplied pixel scales + // (microns/pixel), which may depend on the chip + + float pixelScaleUse, pixelScale1, pixelScale2, pixelScale; + + char *option = psMetadataLookupStr (NULL, recipe, "PSASTRO.COMMON.SCALE.OPTION"); + if (option == NULL) { + psError(PSASTRO_ERR_DATA, false, "no choice set for common scale option\n"); + return false; + } + + bool useExternal = true; + + // find the min or max scale chip + if (!strcasecmp (option, "MIN") || !strcasecmp (option, "MAX")) { + + bool useMax = !strcasecmp (option, "MAX"); + pixelScaleUse = (useMax) ? FLT_MIN : FLT_MAX; + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + + pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]); + pixelScale2 = hypot (chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]); + pixelScale = 0.5*(pixelScale1 + pixelScale2); + + pixelScaleUse = (useMax) ? PS_MAX (pixelScale, pixelScaleUse) : PS_MIN (pixelScale, pixelScaleUse); + } + useExternal = false; + } + + // rescale each chip by the reference scale + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + + psPlaneTransform *toFPA = chip->toFPA; + psPlaneTransform *fromFPA = chip->fromFPA; + + pixelScale1 = hypot (toFPA->x->coeff[1][0], toFPA->x->coeff[0][1]); + pixelScale2 = hypot (toFPA->y->coeff[1][0], toFPA->y->coeff[0][1]); + + if (useExternal) { + pixelScaleUse = getChipPixelScale (chip); + } + + for (int i = 0; i <= toFPA->x->nX; i++) { + for (int j = 0; j <= toFPA->x->nX; j++) { + toFPA->x->coeff[i][j] *= pixelScaleUse/pixelScale1; + toFPA->y->coeff[i][j] *= pixelScaleUse/pixelScale2; + fromFPA->x->coeff[i][j] *= pixelScale1/pixelScaleUse; + fromFPA->y->coeff[i][j] *= pixelScale2/pixelScaleUse; + } + } + + } + return true; +} bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) { @@ -13,5 +85,5 @@ psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); - psPlaneDistortApply (raw->TP, fpa->toTPA, raw->FP, 0.0, 0.0); + psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP); psDeproject (raw->sky, raw->TP, fpa->toSky); }