Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitChip.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitChip.c (revision 12218) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitChip.c (revision 12219) @@ -4,12 +4,14 @@ 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... */ } Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c (revision 12218) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c (revision 12219) @@ -32,4 +32,6 @@ 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,8 +132,12 @@ 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,14 +149,11 @@ 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,123 +164,31 @@ 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); } Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c (revision 12218) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c (revision 12219) @@ -37,4 +37,115 @@ } + +// 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) {