Changeset 21159

Timestamp:

Jan 23, 2009, 3:08:17 PM (17 years ago)

Author:

bills

Message:

add pmAstromLinearizeTransforms to compute a linear fit to the transforms

Location:

trunk/psModules/src/astrom

Files:

• pmAstrometryWCS.c (modified) (2 diffs)
• pmAstrometryWCS.h (modified) (2 diffs)

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -7 +7 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-12-17 19:49:59 $
+ *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-01-24 01:08:17 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -782 +782 @@
 }
 
+static psPlaneTransform *
+linearFitToTransform(psPlaneTransform *trans, psRegion *bounds)
+{
+    int     nSamples = 10;  // 10 samples in each dimension
+
+    double   deltaX = (bounds->x1 - bounds->x0);
+    double   deltaY = (bounds->y1 - bounds->y0);
+
+    psArray *src = psArrayAlloc(nSamples * nSamples);
+    psArray *dst = psArrayAlloc(nSamples * nSamples);
+
+    int k=0;
+    for (int j=0; j<nSamples; j++) {
+        double y = j * deltaY / nSamples;
+        for (int i=0; i<nSamples; i++) {
+            psPlane *s = psPlaneAlloc();
+            s->x = i * deltaX / nSamples;
+            s->y = y;
+            psArraySet(src, k, s);
+            psPlane *d = psPlaneTransformApply(NULL, trans, s);
+            psArraySet(dst, k, d);
+            ++k;
+        }
+    }
+
+    psPlaneTransform *newTrans = psPlaneTransformAlloc(1, 1);
+
+    if (!psPlaneTransformFit(newTrans, src, dst, 0, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
+        return NULL;
+    }
+
+#ifdef COMPARE_TRANS
+    // compare the computed coordintes from this transform with the original
+    psPlane *new = psPlaneAlloc();
+    for (int i=0; i<psArrayLength(dst); i++) {
+        psPlane *d = (psPlane *) psArrayGet(dst, i);
+        psPlane *s = (psPlane *) psArrayGet(src, i);
+
+        new = psPlaneTransformApply(new, newTrans, s);
+
+        printf("%4d %f %f\n", i, 100.*(new->x - d->x)/d->x, 100.*(new->y - d->y)/d->y);
+    }
+    psFree(new);
+#endif
+    psFree(src);
+    psFree(dst);
+
+    return newTrans;
+}
+
+bool pmAstromLinearizeTransforms(pmFPA *fpa, pmChip *chip)
+{
+    psRegion    *chipBounds = pmChipPixels(chip);
+
+    psPlaneTransform *newToFPA = linearFitToTransform(chip->toFPA, chipBounds);
+    if (!newToFPA) {
+        psFree(chipBounds);
+        psError(PS_ERR_UNKNOWN, false, "linear fit for toFPA failed");
+        return false;
+    }
+
+    psFree(chip->toFPA);
+    chip->toFPA = newToFPA;
+
+    psFree(chip->fromFPA);
+    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *chipBounds, 50);
+    if (!chip->fromFPA) {
+        psError(PS_ERR_UNKNOWN, false, "failed to invert linear fit for toFPA");
+        return false;
+    }
+
+    psPlane *chip0 = psPlaneAlloc();
+    chip0->x = 0;
+    chip0->y = 0;
+    psPlane *chip1 = psPlaneAlloc();
+    chip1->x = chipBounds->x1;
+    chip1->y = chipBounds->y1;
+
+    // compute bounding region for fpa
+    psPlane *fpa0 = psPlaneTransformApply(NULL, newToFPA, chip0);
+    psPlane *fpa1 = psPlaneTransformApply(NULL, newToFPA, chip1);
+
+    psRegion *fpaBounds = psRegionAlloc(fpa0->x, fpa1->x, fpa0->y, fpa1->y);
+    psFree(chip0);
+    psFree(chip1);
+    psFree(fpa0);
+    psFree(fpa1);
+
+    psPlaneTransform *newToTPA = linearFitToTransform(fpa->toTPA, fpaBounds);
+    if (!newToTPA) {
+        psError(PS_ERR_UNKNOWN, false, "failed to perform linear fit to toTPA");
+        psFree(fpaBounds);
+        return false;
+    }
+    psFree(fpa->toTPA);
+    fpa->toTPA = newToTPA;
+
+    // XXX: is this region ok?
+    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *fpaBounds, 50);
+    if (!fpa->fromTPA) {
+        psError(PS_ERR_UNKNOWN, false, "failed to invert linear fit to toTPA");
+        return false;
+    }
+
+    fpa->toSky->type = PS_PROJ_TAN;
+
+    psFree(chipBounds);
+    psFree(fpaBounds);
+
+    return true;
+}
+
 static void pmAstromWCSFree (pmAstromWCS *wcs)
 {

trunk/psModules/src/astrom/pmAstrometryWCS.h

@@ -4 +4 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-03-21 22:00:49 $
+ * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-01-24 01:08:17 $
  * Copyright 2006 Institute for Astronomy, University of Hawaii
  */
@@ -62 +62 @@
 bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, double tol);
 
+bool pmAstromLinearizeTransforms(pmFPA *fpa, pmChip *);
+
 // move to pslib
 psPlaneDistort *psPlaneDistortIdentity (int order);