Index: /branches/czw_branch/20160809/psastro/src/psastroAstromGuess.c =================================================================== --- /branches/czw_branch/20160809/psastro/src/psastroAstromGuess.c (revision 39685) +++ /branches/czw_branch/20160809/psastro/src/psastroAstromGuess.c (revision 39686) @@ -75,5 +75,6 @@ pmFPA *fpa = input->fpa; - if (DEBUG) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa); + // this call only works if we have loaded a model : seg fault if not + if (DEBUG && useModel) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa); // load mosaic-level astrometry? @@ -91,4 +92,10 @@ if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; } + + if (!chip->toFPA || !chip->fromFPA) { + char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + fprintf (stderr, "no astrom model for %s, skipping\n", name); + continue; + } if (newFPA) { @@ -312,4 +319,10 @@ while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } + + if (!chip->toFPA || !chip->fromFPA) { + char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + fprintf (stderr, "no astrom model for %s, skipping\n", name); + continue; + } // XXX we are currently inconsistent with marking the good vs the bad data Index: /branches/czw_branch/20160809/psastro/src/psastroDemoDump.c =================================================================== --- /branches/czw_branch/20160809/psastro/src/psastroDemoDump.c (revision 39685) +++ /branches/czw_branch/20160809/psastro/src/psastroDemoDump.c (revision 39686) @@ -234,4 +234,10 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } + if (!chip->toFPA || !chip->fromFPA) { + char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + fprintf (stderr, "no astrom model for %s, skipping\n", name); + continue; + } + // XXX write out the four corners for a test psRegion *region = pmChipPixels (chip); Index: /branches/czw_branch/20160809/psastro/src/psastroModelAdjust.c =================================================================== --- /branches/czw_branch/20160809/psastro/src/psastroModelAdjust.c (revision 39685) +++ /branches/czw_branch/20160809/psastro/src/psastroModelAdjust.c (revision 39686) @@ -52,4 +52,9 @@ return false; } + float refChipAngleNominal = PS_RAD_DEG*psMetadataLookupF32 (&status, recipe, "PSASTRO.MODEL.REF.CHIP.ANGLE"); + if (!refChipName) { + psError(PS_ERR_IO, true, "reference chip is missing from recipe"); + return false; + } // get reference chip from name @@ -87,5 +92,11 @@ // get the current posangle of the ref chip float chipAngle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]); - fprintf (stderr, "chipAngle: %f\n", chipAngle*PS_DEG_RAD); + + // chipAngle should be refChipAngleNominal @ POSANGLE = 0.0 + float posAngleOffset = chipAngle - refChipAngleNominal; + + fprintf (stderr, "chipAngle is: %f, at PA = 0.0, it should be %f, rotating model by %f\n", + chipAngle*PS_DEG_RAD, refChipAngleNominal*PS_DEG_RAD, posAngleOffset*PS_DEG_RAD); + // psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle); @@ -99,5 +110,5 @@ psRegion *region = pmChipPixels (chip); - psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, chipAngle); + psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, posAngleOffset); psFree (chip->toFPA); chip->toFPA = toFPA; Index: /branches/czw_branch/20160809/psastro/src/psastroMosaicAstrom.c =================================================================== --- /branches/czw_branch/20160809/psastro/src/psastroMosaicAstrom.c (revision 39685) +++ /branches/czw_branch/20160809/psastro/src/psastroMosaicAstrom.c (revision 39686) @@ -15,4 +15,5 @@ bool psastroMosaicFit (pmFPA *fpa, psMetadata *recipe, const char *rootname, int pass); +bool psastroProjectionRefit (pmFPA *fpa, psMetadata *recipe); // XXX require this fpa to have multiple chip extensions and a PHU? @@ -42,4 +43,9 @@ int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); if (!status) psAbort ("missing config value"); + + // XXX if projection is not TAN, fit the measured projection here and combine toTPA/fromTPA functions + if ((fpa->toSky->type != PS_PROJ_TAN) && (fpa->toSky->type != PS_PROJ_DIS)) { + if (!psastroProjectionRefit (fpa, recipe)) psAbort ("failed to refit distortion"); + } // this should be in a loop with nIter = @@ -154,2 +160,134 @@ return true; } + +// we have a fpa->toSky projection which is NOT of type TAN along with a toTPA which is +// the Identity transform. we want to convert this to TAN projection plus a non-identity +// transformation to take the non-TAN components. + +// we are going to generate a set of FP->(x,y) using the current projection + fromTPA +// transformations along with a set of TP->(x,y) values using the desired TAN projection, +// then we will fit TP vs FP + +bool psastroProjectionRefit (pmFPA *fpa, psMetadata *recipe) { + + bool status; + + // 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 mosaic distortion fit order\n"); + return false; + } + + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + pmFPAview *view = pmFPAviewAlloc (0); + + psVector *L = psVectorAllocEmpty (1000, PS_TYPE_F32); + psVector *M = psVectorAllocEmpty (1000, PS_TYPE_F32); + + psVector *P = psVectorAllocEmpty (1000, PS_TYPE_F32); + psVector *Q = psVectorAllocEmpty (1000, PS_TYPE_F32); + + psProjection *toSkyTan = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, fpa->toSky->Xs, fpa->toSky->Ys, PS_PROJ_TAN); + + float xMin = NAN; + float xMax = NAN; + float yMin = NAN; + float yMax = NAN; + + bool firstObject = true; + + // 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 *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); + if (refstars == NULL) continue; + psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); + + psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); + if (matches == NULL) continue; + + // we are looking over the matched refstars only to be sure we are + // covering the valid space of the projection + transformation, and not + // beyond + + for (int i = 0; i < matches->n; i++) { + pmAstromMatch *match = matches->data[i]; + pmAstromObj *ref = refstars->data[match->ref]; + + psPlane fpOld, tpNew, tpOld; + + psProject (&tpOld, ref->sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection + psPlaneTransformApply (&fpOld, fpa->fromTPA, &tpOld); + + psProject (&tpNew, ref->sky, toSkyTan); // find the focal-plane coord of this RA,DEC coord using the ref chip projection + + psVectorAppend (L, fpOld.x); + psVectorAppend (M, fpOld.y); + + psVectorAppend (P, tpNew.x); + psVectorAppend (Q, tpNew.y); + + if (firstObject) { + xMin = xMax = fpOld.x; + yMin = yMax = fpOld.y; + firstObject = false; + } + + xMin = PS_MIN (xMin, fpOld.x); + xMax = PS_MAX (xMax, fpOld.x); + yMin = PS_MIN (yMin, fpOld.y); + yMax = PS_MAX (yMax, fpOld.y); + } + } + } + } + + // the original transforms are (1, 1), but we will need higher order to fit the distortions + psFree (fpa->toTPA); + + // the original transforms are (1, 1), but we will need higher order to fit the distortions + 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->coeffMask[i][j] = PS_POLY_MASK_SET; + fpa->toTPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; + } + } + } + + psVectorFitPolynomial2D (fpa->toTPA->x, NULL, 0, P, NULL, L, M); + psVectorFitPolynomial2D (fpa->toTPA->y, NULL, 0, Q, NULL, L, M); + + psRegion fitRegion = psRegionSet (xMin, xMax, yMin, yMax); + + // psPlaneTransformInvert will generate a new fromTPA with order to match toTPA + fpa->fromTPA = psPlaneTransformInvert (fpa->fromTPA, fpa->toTPA, fitRegion, 100); + + psFree (fpa->toSky); + fpa->toSky = toSkyTan; + + psFree (L); + psFree (M); + psFree (P); + psFree (Q); + + psFree (view); + return true; +} + Index: /branches/czw_branch/20160809/psastro/src/psastroMosaicGradients.c =================================================================== --- /branches/czw_branch/20160809/psastro/src/psastroMosaicGradients.c (revision 39685) +++ /branches/czw_branch/20160809/psastro/src/psastroMosaicGradients.c (revision 39686) @@ -79,5 +79,5 @@ int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER"); if (!status) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n"); + psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic distortion fit order\n"); psFree (gradients); psFree (view); Index: /branches/czw_branch/20160809/psastro/src/psastroMosaicOneChip.c =================================================================== --- /branches/czw_branch/20160809/psastro/src/psastroMosaicOneChip.c (revision 39685) +++ /branches/czw_branch/20160809/psastro/src/psastroMosaicOneChip.c (revision 39686) @@ -100,19 +100,19 @@ } - // XXX allow statitic to be set by the user + // XXX allow statistic to be set by the user // only clip if we are fitting the chip parameters. psStats *fitStats = NULL; -// if (order == 0) { -// 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"); -// } else { - 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"); -// } + if (FALSE && (order == 0)) { + 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"); + } else { + 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); + pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats, recipe); if (!results) { psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n"); Index: /branches/czw_branch/20160809/psastro/src/psastroOneChipFit.c =================================================================== --- /branches/czw_branch/20160809/psastro/src/psastroOneChipFit.c (revision 39685) +++ /branches/czw_branch/20160809/psastro/src/psastroOneChipFit.c (revision 39686) @@ -124,5 +124,5 @@ // improved fit for astrometric terms - results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); + results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats, recipe); if (!results) { psLogMsg ("psastro", 3, "failed to perform the matched fit\n");