Index: /trunk/pswarp/src/pswarp.h =================================================================== --- /trunk/pswarp/src/pswarp.h (revision 36834) +++ /trunk/pswarp/src/pswarp.h (revision 36835) @@ -30,7 +30,10 @@ #define PSASTRO_RECIPE "PSASTRO" ///< Name of the recipe to use -#define PSWARP_ANALYSIS_GOODPIX "PSWARP.GOODPIX" ///< Name for number of good pixels in analysis metadata +#define PSWARP_ANALYSIS_GOODPIX "PSWARP.GOODPIX" ///< Name for number of good pixels in analysis metadata #define PSWARP_ANALYSIS_COVARIANCES "PSWARP.COVARIANCES" ///< Name for covariance matrices on analysis MD -#define PSWARP_ANALYSIS_JACOBIAN "PSWARP.JACOBIAN" ///< Name for Jacobian on analysis MD +#define PSWARP_ANALYSIS_JACOBIAN "PSWARP.JACOBIAN" ///< Name for Jacobian on analysis MD +#define PSWARP_ANALYSIS_BACKMAPS "PSWARP.BACKMAPS" ///< Name for array of backwarps maps on analysis MD +#define PSWARP_ANALYSIS_CHIPNAMES "PSWARP.CHIPNAMES" ///< Name for array of input chipnames on analysis MD +#define PSWARP_ANALYSIS_CHIPREGIONS "PSWARP.CHIPREGIONS" ///< Name for array of input chipnames on analysis MD /** Index: /trunk/pswarp/src/pswarpLoop.c =================================================================== --- /trunk/pswarp/src/pswarpLoop.c (revision 36834) +++ /trunk/pswarp/src/pswarpLoop.c (revision 36835) @@ -156,6 +156,6 @@ if (!pswarpMakePSF (config, output, stats)) { - psError(psErrorCodeLast(), false, "problem generating PSF."); - goto FAIL; + psError(psErrorCodeLast(), false, "problem generating PSF."); + goto FAIL; } @@ -163,8 +163,10 @@ return true; - FAIL: +FAIL: psFree (view); return false; } + +bool pswarpGetBackTransform (pmReadout *tgt, pmReadout *src); // once the output fpa elements have been built, loop over the fpa and generate stats @@ -188,4 +190,7 @@ while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) { pswarpTransformReadout (readout, input, config, backgroundWarp); + + // determine a (rough) linear WCS from readout to input + pswarpGetBackTransform (readout, input); } } @@ -194,2 +199,118 @@ } +bool pswarpGetBackTransform (pmReadout *tgt, pmReadout *src) { + + bool status; + + // tgt is the pswarp output target image, src is the input source image + // here we are generating a linear transformation from tgt -> src + + // the map is defined for coordinates in the image parent frame. + psRegion imageRegion = psRegionSet (0,0,0,0); + imageRegion = psRegionForImage (src->image, imageRegion); + psString imageRegionString = psRegionToString (imageRegion); + + // generate a set of points (X,Y)_tgt = (X,Y)_src and do the 2D fit? + // save a 2D polynomial, not WCS? + + // use a map with only linear terms in x,y + psPlaneTransform *map = psPlaneTransformAlloc(1,1); + map->x->coeffMask[1][1] = PS_POLY_MASK_BOTH; + map->y->coeffMask[1][1] = PS_POLY_MASK_BOTH; + + pmCell *cell = NULL; + + cell = src->parent; + pmChip *chipSrc = cell->parent; + pmFPA *fpaSrc = chipSrc->parent; + + cell = tgt->parent; + pmChip *chipTgt = cell->parent; + pmFPA *fpaTgt = chipTgt->parent; + + psPlane *srcPos = psPlaneAlloc(); + psPlane *tgtPos = psPlaneAlloc(); + + psPlane *FP = psPlaneAlloc(); + psPlane *TP = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + + psVector *tgtX = psVectorAllocEmpty (121, PS_TYPE_F32); + psVector *tgtY = psVectorAllocEmpty (121, PS_TYPE_F32); + psVector *srcX = psVectorAllocEmpty (121, PS_TYPE_F32); + psVector *srcY = psVectorAllocEmpty (121, PS_TYPE_F32); + + int dX = src->image->numCols / 10.0; + int dY = src->image->numRows / 10.0; + for (int ix = 0; ix < src->image->numCols; ix += dX) { + for (int iy = 0; iy < src->image->numRows; iy += dY) { + + srcPos->x = ix; + srcPos->y = iy; + + psPlaneTransformApply(FP, chipSrc->toFPA, srcPos); + psPlaneTransformApply (TP, fpaSrc->toTPA, FP); + psDeproject (sky, TP, fpaSrc->toSky); + + psProject (TP, sky, fpaTgt->toSky); + psPlaneTransformApply (FP, fpaTgt->fromTPA, TP); + psPlaneTransformApply (tgtPos, chipTgt->fromFPA, FP); + + psVectorAppend (srcX, srcPos->x); + psVectorAppend (srcY, srcPos->y); + psVectorAppend (tgtX, tgtPos->x); + psVectorAppend (tgtY, tgtPos->y); + } + } + + psVectorFitPolynomial2D(map->x, NULL, 0, srcX, NULL, tgtX, tgtY); + psVectorFitPolynomial2D(map->y, NULL, 0, srcY, NULL, tgtX, tgtY); + + // for each output image, we want to add headers corresponding to all input images + // save an array of the transformations on the tgt analysis and an array of the input chip names + psArray *backmaps = psMetadataLookupPtr (&status, tgt->analysis, PSWARP_ANALYSIS_BACKMAPS); + if (!backmaps) { + backmaps = psArrayAllocEmpty (4); + psMetadataAddArray (tgt->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_BACKMAPS, 0, "backwards maps", backmaps); + psFree (backmaps); // can free here since we put a copy on tgt->analysis + } + psArrayAdd (backmaps, 4, map); + + // we also need to save the corresponding chip name for this map + char *chipname = psMetadataLookupStr (&status, chipSrc->concepts, "CHIP.NAME"); + psArray *chipnames = psMetadataLookupPtr (&status, tgt->analysis, PSWARP_ANALYSIS_CHIPNAMES); + if (!chipnames) { + chipnames = psArrayAllocEmpty (4); + psMetadataAddArray (tgt->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_CHIPNAMES, 0, "source images", chipnames); + psFree (chipnames); // can free here since we put a copy on tgt->analysis + } + psArrayAdd (chipnames, 4, chipname); + + // we also need to save the corresponding chip name for this map + psArray *chipRegions = psMetadataLookupPtr (&status, tgt->analysis, PSWARP_ANALYSIS_CHIPREGIONS); + if (!chipRegions) { + chipRegions = psArrayAllocEmpty (4); + psMetadataAddArray (tgt->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_CHIPREGIONS, 0, "source images", chipRegions); + psFree (chipRegions); // can free here since we put a copy on tgt->analysis + } + psArrayAdd (chipRegions, 4, imageRegionString); + + fprintf (stderr, "warp to %s X = %f + %f x + %f y\n", chipname, map->x->coeff[0][0], map->x->coeff[1][0], map->x->coeff[0][1]); + fprintf (stderr, "warp to %s Y = %f + %f x + %f y\n", chipname, map->y->coeff[0][0], map->y->coeff[1][0], map->y->coeff[0][1]); + + psFree (map); + psFree (tgtX); + psFree (tgtY); + psFree (srcX); + psFree (srcY); + + psFree (srcPos); + psFree (tgtPos); + psFree (FP); + psFree (TP); + psFree (sky); + + psFree (imageRegionString); + + return true; +} Index: /trunk/pswarp/src/pswarpUpdateMetadata.c =================================================================== --- /trunk/pswarp/src/pswarpUpdateMetadata.c (revision 36834) +++ /trunk/pswarp/src/pswarpUpdateMetadata.c (revision 36835) @@ -78,4 +78,40 @@ psFree(md5string); psFree(headerName); + + bool status; + psString keyword = NULL; + psString mapstring = NULL; + + // we (should) have an array of chipnames on the analysis + psArray *chipnames = psMetadataLookupPtr(&status, readout->analysis, PSWARP_ANALYSIS_CHIPNAMES); + for (int i = 0; chipnames && (i < chipnames->n); i++) { + psStringAppend (&keyword, "SRC_%04d", i); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", chipnames->data[i]); + psFree (keyword); + } + + // we (should) have an array of chipnames on the analysis + psArray *chipRegions = psMetadataLookupPtr(&status, readout->analysis, PSWARP_ANALYSIS_CHIPREGIONS); + for (int i = 0; chipRegions && (i < chipRegions->n); i++) { + psStringAppend (&keyword, "SEC_%04d", i); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", chipRegions->data[i]); + psFree (keyword); + } + + // we (should) also have an array of backwards maps (warp->chip) for the above chips + psArray *backmaps = psMetadataLookupPtr(&status, readout->analysis, PSWARP_ANALYSIS_BACKMAPS); // + for (int i = 0; backmaps && (i < backmaps->n); i++) { + psPlaneTransform *map = backmaps->data[i]; + psStringAppend (&keyword, "MPX_%04d", i); + psStringAppend (&mapstring, "[%8.2f,%8.4f,%8.4f]", map->x->coeff[0][0], map->x->coeff[1][0], map->x->coeff[0][1]); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "warp to image map", mapstring); + psFree (keyword); + psFree (mapstring); + psStringAppend (&keyword, "MPY_%04d", i); + psStringAppend (&mapstring, "[%8.2f,%8.4f,%8.4f]", map->y->coeff[0][0], map->y->coeff[1][0], map->y->coeff[0][1]); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "warp to image map", mapstring); + psFree (keyword); + psFree (mapstring); + } }