Changeset 36835 for trunk/pswarp/src/pswarpLoop.c

Timestamp:

Jun 7, 2014, 6:38:38 AM (12 years ago)

Author:

eugene

Message:

add pswarp backwards transformations so we can map warp pixels to chip pixels

File:

• trunk/pswarp/src/pswarpLoop.c (modified) (4 diffs)

trunk/pswarp/src/pswarpLoop.c

@@ -156 +156 @@
 
     if (!pswarpMakePSF (config, output, stats)) {
-      psError(psErrorCodeLast(), false, "problem generating PSF.");
-      goto FAIL;
+        psError(psErrorCodeLast(), false, "problem generating PSF.");
+        goto FAIL;
     }
 
@@ -163 +163 @@
     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 +190 @@
             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 +199 @@
 }
 
+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;
+}