Changeset 12219

Timestamp:

Mar 4, 2007, 11:56:49 AM (19 years ago)

Author:

eugene

Message:

fixing fits

Location:

branches/dvo-mods-2007-02/Ohana/src/relastro/src

Files:

• FitChip.c (modified) (1 diff)
• fitpoly.c (modified) (4 diffs)
• mkpolyterm.c (modified) (1 diff)

branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitChip.c

@@ -4 +4 @@
 
   int i;
+  CoordFit *fit;
 
-  fit_init (coords[0].Npolyterms);
+  fit = fit_init (coords[0].Npolyterms);
   for (i = 0; i < Nmatch; i++) {
     if (raw[i].mask) continue;
-    fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M);
+    fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0);
   }
-  fit_eval ();
-  fit_apply_coords (coords);
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
   /* FitChip and FitSimple update the coords in different ways? maybe not... */
 }

branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c

@@ -32 +32 @@
   ALLOCATE (fit[0].ysum, double *, fit[0].Nterms);
   for (i = 0; i < fit[0].Nterms; i++) {
+    ALLOCATE (fit[0].xsum[i], double, fit[0].Nterms);
+    ALLOCATE (fit[0].ysum[i], double, fit[0].Nterms);
     bzero (fit[0].xsum[i], fit[0].Nterms*sizeof(double));
     bzero (fit[0].ysum[i], fit[0].Nterms*sizeof(double));
@@ -130 +132 @@
   modfit = CoordsSetCenter (fit, xo, yo);
 
-  // set pc1_1, pc1_2, pc2_1, pc2_2
+  // set crpix1,2
+  fit_to_coordterms (coords, modfit, 0, 0);
+
+  // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0)
   fit_to_coordterms (coords, modfit, 0, 1);
   fit_to_coordterms (coords, modfit, 1, 0);
 
+  // set the polyterm elements 
   fit_to_coordterms (coords, modfit, 0, 2);
   fit_to_coordterms (coords, modfit, 1, 1);
@@ -143 +149 @@
   fit_to_coordterms (coords, modfit, 3, 0);
 
-  /* get the correct vector entries for the linear terms */
-  coords[0].crval1 = vector[N][0];  
-  coords[0].crval2 = vector[N][1];
+  /* we do not modify crval1,2: these are kept at the default values */
 
-  N = mkvector (1, 0, NORDER);
-  c11 = vector[N][0];  
-  c21 = vector[N][1];
-  N = mkvector (0, 1, NORDER);
-  c12 = vector[N][0];  
-  c22 = vector[N][1];
+  /* normalize pc11,etc */
+  c11 = coords[0].pc1_1;
+  c21 = coords[0].pc1_1;
+  c12 = coords[0].pc1_1;
+  c22 = coords[0].pc1_1;
   coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21));
   R = 1 / coords[0].cdelt1;
@@ -161 +164 @@
   coords[0].pc2_2  = c22*R;
 
-  coords[0].crpix1 = 0;
-  coords[0].crpix2 = 0;
-
-  coords[0].Npolyterms = NORDER;
-  strcpy (coords[0].ctype, "DEC--PLY");
+  /* keep the order and type from initial values */
 }
 
-/*
-  if we have just linear terms, the following holds for crpix1,2:
-  D = R / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
-  coords[0].crpix1 = D * (coords[0].pc1_2*c20 - coords[0].pc2_2*c10);
-  coords[0].crpix2 = D * (coords[0].pc2_1*c10 - coords[0].pc1_1*c20);
-*/
-
-/* NORDER is order of gradient fit : Npolyterms is Norder + 1 */
-void fit_apply_grads (Coords *distort, Coords *project, int term) {
-
-  int i, j, Np, Nv1, Nv2;
-
-  distort[0].Npolyterms = NORDER + 1;
-  if (distort[0].Npolyterms > 1) strcpy (distort[0].ctype, "DEC--PLY");
-
-  for (i = 0; i < NPOWER + 1; i++) {
-    for (j = 0; j < (NPOWER + 1 - i); j++) {
-      
-      if (i + j < 2) continue;
-      Np  = mkpolyterm (i, j);
-      Nv1 = mkvector (i-1, j, NORDER);
-      Nv2 = mkvector (i, j-1, NORDER);
-
-      /** why do we have the negative sign? **/
-      if (j == 0) {
-        distort[0].polyterms[Np][term] = vector[Nv1][0] / i;
-      }
-      if (i == 0) {
-        distort[0].polyterms[Np][term] = vector[Nv2][1] / j;
-      }
-      if ((i > 0) && (j > 0)) {
-        distort[0].polyterms[Np][term] = 0.5*(vector[Nv1][0] / i + vector[Nv2][1] / j);
-      }
-    }
-  }
-
-  Nv1 = mkvector (0, 0, NORDER);
-  if (term == 0) {
-    project[0].pc1_1 = project[0].pc1_1 * (1 + vector[Nv1][0]);
-    project[0].pc1_2 = project[0].pc1_2 * (1 + vector[Nv1][1]);
-  } else {
-    project[0].pc2_1 = project[0].pc2_1 * (1 + vector[Nv1][0]);
-    project[0].pc2_2 = project[0].pc2_2 * (1 + vector[Nv1][1]);
-  }
-}
-
-void fit_correct_grads (Gradients *in, Gradients *out, int term) {
-
-  int i, k, m, n;
-  double x, y, dx, dy, dz1, dz2;
-
-  for (i = 0; i < in[0].Npts; i++) {
-    
-    dx = in[0].Lo[i];
-    dy = in[0].Mo[i];
-    dz1 = dz2 = 0.0;
-
-    k = 0;
-    x = y = 1;
-    for (m = 0; m < NPOWER; m++) {
-      x = y;
-      for (n = 0; n < NPOWER - m; n++, k++) {
-        dz1 += vector[k][0]*x;
-        dz2 += vector[k][1]*x;
-        x = x * dx / SCALE;
-      }
-      y = y * dy / SCALE;
-    }
-
-    out[0].Lo[i] = dx;
-    out[0].Mo[i] = dy;
-    if (term == 0) {
-      out[0].dPdL[i] = in[0].dPdL[i] - dz1;
-      out[0].dPdM[i] = in[0].dPdM[i] - dz2;
-    } else {
-      out[0].dQdL[i] = in[0].dQdL[i] - dz1;
-      out[0].dQdM[i] = in[0].dQdM[i] - dz2;
-    }
-  }
-}
-
-int mkvector (int n, int m, int norder) {
-  
-  int i, N;
-  
-  N = 0;
-  for (i = 0; i < m; i++) {
-    N += (norder - i + 1);
-  }
-  N += n;
-  return (N);
-}
-
-void fit_free () {
+void fit_free (CoordFit *fit) {
 
   int i;
 
-  for (i = 0; i < NTERM; i++) {
-    free (sum[i]);
-    free (xsum[i]);
-    free (ysum[i]);
+  for (i = 0; i < fit[0].Nterms; i++) {
+    free (fit[0].xsum[i]);
+    free (fit[0].ysum[i]);
+  }      
+  free (fit[0].xsum);
+  free (fit[0].ysum);
+
+  for (i = 0; i < fit[0].Nsums; i++) {
+    free (fit[0].sum[i]);
   }
-  free (sum);
-  free (xsum);
-  free (ysum);
+  free (fit[0].sum);
 
-  for (i = 0; i < NPARS; i++) {
-    free (matrix[i]);
-    free (vector[i]);
+  for (i = 0; i < fit[0].Nelems; i++) {
+    free (fit[0].matrix[i]);
+    free (fit[0].vector[i]);
   }
-  free (matrix);
-  free (vector);
+  free (fit[0].matrix);
+  free (fit[0].vector);
+
+  free (fit);
 }
   

branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c

@@ -37 +37 @@
 }
 
+
+// XXXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ?
+
+CoordsGetCenter () {
+
+// given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y
+// for which L,M = 0,0. tol is the allowed error on x,y.
+psPlane *psPlaneTransformGetCenter (psPlaneTransform *trans, double tol)
+{
+
+    // crpix1,2 = X,Y(crval1,2)
+    // start with linear solution for Xo,Yo:
+    double R  = (trans->x->coeff[1][0]*trans->y->coeff[0][1] - trans->x->coeff[0][1]*trans->y->coeff[1][0]);
+    double Xo = (trans->y->coeff[0][0]*trans->x->coeff[0][1] - trans->x->coeff[0][0]*trans->y->coeff[0][1])/R;
+    double Yo = (trans->x->coeff[0][0]*trans->y->coeff[1][0] - trans->y->coeff[0][0]*trans->x->coeff[1][0])/R;
+
+    // iterate to actual solution: requires small non-linear terms
+    if (trans->x->nX > 1) {
+        psPolynomial2D *XdX = psPolynomial2D_dX(NULL, trans->x);
+        psPolynomial2D *XdY = psPolynomial2D_dY(NULL, trans->x);
+
+        psPolynomial2D *YdX = psPolynomial2D_dX(NULL, trans->y);
+        psPolynomial2D *YdY = psPolynomial2D_dY(NULL, trans->y);
+
+        psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32);
+        psVector *Beta = psVectorAlloc (2, PS_DATA_F32);
+
+        /* this loop uses the Newton-Raphson method to solve for Xo,Yo
+        * it needs the high order terms to be small 
+        * Xo,Yo are in pixels;
+        */
+        double dPos = tol + 1;
+        for (int i = 0; (dPos > tol) && (i < 20); i++) {
+            // NOTE: order for Alpha is: [y][x]
+            Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo);
+            Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
+            Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
+            Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
+
+            Beta->data.F32[0] = psPolynomial2DEval (trans->x, Xo, Yo);
+            Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo);
+
+            psMatrixGJSolveF32 (Alpha, Beta);
+
+            Xo -= Beta->data.F32[0];
+            Yo -= Beta->data.F32[1];
+            dPos = hypot(Beta->data.F32[0], Beta->data.F32[1]);
+        }
+        psFree (Alpha);
+        psFree (Beta);
+        psFree (XdX);
+        psFree (XdY);
+        psFree (YdX);
+        psFree (YdY);
+    }
+    psPlane *center = psPlaneAlloc ();
+    center->x = Xo;
+    center->y = Yo;
+
+    return center;
+}
+
+// convert a transformation L(x,y) to L'(x-xo,y-yo)
+psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, double Xo, double Yo)
+{
+
+    // validate fit order
+
+    if (output == NULL) {
+        output = psPlaneTransformAlloc(input->x->nX, input->x->nY);
+    }
+
+    /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
+     * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and 
+     * noting that the constant term in each is the coefficient in the case of L(x,y) and is the 
+     * value of L'(-xo,-yo) in the second case.
+     */
+
+    psPolynomial2D *tmp;
+
+    psPolynomial2D *xPx = psPolynomial2DCopy (NULL, input->x);
+    psPolynomial2D *yPx = psPolynomial2DCopy (NULL, input->y);
+
+    for (int i = 0; i <= input->x->nX; i++) {
+        psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
+        psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
+        for (int j = 0; j <= input->x->nY; j++) {
+            output->x->mask[i][j] = input->x->mask[i][j];
+            output->y->mask[i][j] = input->y->mask[i][j];
+            output->x->coeff[i][j] = (output->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
+            output->y->coeff[i][j] = (output->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
+
+            // take the next derivative wrt y, catch output (is NULL on last pass)
+            tmp = psPolynomial2D_dY(NULL, xPy);
+            psFree (xPy);
+            xPy = tmp;
+            tmp = psPolynomial2D_dY(NULL, yPy);
+            psFree (yPy);
+            yPy = tmp;
+        }
+        // take the next derivative wrt x, catch output (is NULL on last pass)
+        tmp = psPolynomial2D_dX(NULL, xPx);
+        psFree (xPx);
+        xPx = tmp;
+        tmp = psPolynomial2D_dX(NULL, yPx);
+        psFree (yPx);
+        yPx = tmp;
+    }
+    return output;
+}
+
 int mkpolyterm (int n, int m) {