Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 6387) +++ /trunk/psLib/src/astro/psCoord.c (revision 6388) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.109 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-02-08 21:52:05 $ +* @version $Revision: 1.110 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-02-09 00:52:18 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -28,4 +28,5 @@ #include "psErrorText.h" #include "psMatrix.h" +#include "psMinimizePolyFit.h" #include #include @@ -461,6 +462,9 @@ out = psSphereAlloc(); } else { + if ( !psMemCheckType(PS_DATA_SPHERE, outSphere) ) { + psError(PS_ERR_UNKNOWN, true, "outSphere must be of type *psSphere.\n"); + return(NULL); + } out = outSphere; - // XXX: Do a memory checkpointer. } @@ -565,10 +569,4 @@ for (psS32 t2x = 0 ; t2x < (1 + trans2->nX) ; t2x++) { for (psS32 t2y = 0 ; t2y < (1 + trans2->nY) ; t2y++) { - /* Possible debug-only macro which checks these coords? - if ((t1x+t2x) >= orderX) - printf("BAD 1\n"); - if ((t1y+t2y) >= orderY) - printf("BAD 2\n"); - */ out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]); } @@ -625,5 +623,4 @@ // Allocate the new psPlaneTransform, if necessary. // - // XXX: rename, or verify, or recode, after the poly norder/nterm change. psPlaneTransform *myPT = NULL; @@ -639,6 +636,6 @@ // Initialize the new psPlaneTransform, if necessary. // - for (psS32 i = 0 ; i < orderX ; i++) { - for (psS32 j = 0 ; j < orderY ; j++) { + for (psS32 i = 0 ; i < orderX+1 ; i++) { + for (psS32 j = 0 ; j < orderY+1 ; j++) { myPT->x->coeff[i][j] = 0.0; myPT->x->mask[i][j] = 0; @@ -769,6 +766,5 @@ psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip) -XXX: What about nRejIter? Iterations? -XXX: This code has problems with data that corresponds to a non-linear fit. +XXX: This code ignores nRejIter and sigmaClip. *****************************************************************************/ bool psPlaneTransformFit( @@ -783,112 +779,34 @@ PS_ASSERT_PTR_NON_NULL(dest, NULL); - psBool rc = true; + // + // Create the x and y vectors for the psVectorFitPolynomial2D() function. + // psS32 numCoords = PS_MIN(source->n, dest->n); - psS32 order = PS_MAX(trans->x->nX, trans->x->nY); - order = PS_MAX(order, trans->y->nX); - order = PS_MAX(order, trans->y->nY); - // XXX: Verify, or recode, after the poly norder/nterm change. - order++; - - // - // Create fake polynomial to use in evaluation - // - // XXX: Verify this (poly order/nterm change) - psPolynomial2D *fakePoly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order-1, order-1); - for (int i = 0; i < order; i++) { - for (int j = 0; j < order; j++) { - fakePoly->coeff[i][j] = 1.0; - fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate - } - } - - // - // Initialize the matrix and vectors - // - psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients - psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution - psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x - psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y - for (psS32 i = 0; i < nCoeff; i++) { - for (psS32 j = 0; j < nCoeff; j++) { - matrix->data.F64[i][j] = 0.0; - } - xVector->data.F64[i] = 0.0; - yVector->data.F64[i] = 0.0; - } - - // - // Iterate over the grid points - // + psVector *xIn = psVectorAlloc(numCoords, PS_TYPE_F64); + psVector *yIn = psVectorAlloc(numCoords, PS_TYPE_F64); + psVector *xOut = psVectorAlloc(numCoords, PS_TYPE_F64); + psVector *yOut = psVectorAlloc(numCoords, PS_TYPE_F64); for (psS32 g = 0; g < numCoords; g++) { - // Iterate over the polynomial coefficients, accumulating the matrix and vectors - - for (psS32 i = 0, ijIndex = 0; i < order; i++) { - for (psS32 j = 0; j < order - i; j++, ijIndex++) { - fakePoly->mask[i][j] = 0; - psF64 xIn = ((psPlane *) source->data[g])->x; - psF64 yIn = ((psPlane *) source->data[g])->y; - psF64 xOut = ((psPlane *) dest->data[g])->x; - psF64 yOut = ((psPlane *) dest->data[g])->y; - psF64 ijPoly = psPolynomial2DEval(fakePoly, xIn, yIn); - fakePoly->mask[i][j] = 1; - - for (psS32 m = 0, mnIndex = 0; m < order; m++) { - for (psS32 n = 0; n < order - m; n++, mnIndex++) { - fakePoly->mask[m][n] = 0; - psF64 mnPoly = psPolynomial2DEval(fakePoly, xIn, yIn); - fakePoly->mask[m][n] = 1; - - matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly; - } - } - - xVector->data.F64[ijIndex] += ijPoly * xOut; - yVector->data.F64[ijIndex] += ijPoly * yOut; - } - } - } - - // - // Solution via LU Decomposition - // - psVector *xSolution; - psVector *ySolution; - psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition - psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix - if (luMatrix == NULL) { - rc = false; - psError(PS_ERR_UNKNOWN, true, "psMatrixLUD() returned NULL. Returning FALSE.\n"); - } else { - xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x - ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y - - if ((xSolution == NULL) || (ySolution == NULL)) { - rc = false; - psError(PS_ERR_UNKNOWN, true, "psMatrixLUSolve() returned NULL. Returning FALSE.\n"); - } else { - // - // Stuff coefficients into transformation - // - for (psS32 i = 0, ijIndex = 0; i < order; i++) { - for (psS32 j = 0; j < order - i; j++, ijIndex++) { - trans->x->coeff[i][j] = xSolution->data.F64[ijIndex]; - trans->y->coeff[i][j] = ySolution->data.F64[ijIndex]; - } - } - } - } - - psFree(fakePoly); - psFree(permutation); - psFree(luMatrix); - psFree(xSolution); - psFree(ySolution); - psFree(matrix); - psFree(xVector); - psFree(yVector); - - return(rc); -} + xIn->data.F64[g] = ((psPlane *) source->data[g])->x; + yIn->data.F64[g] = ((psPlane *) source->data[g])->y; + xOut->data.F64[g] = ((psPlane *) dest->data[g])->x; + yOut->data.F64[g] = ((psPlane *) dest->data[g])->y; + } + + trans->x = psVectorFitPolynomial2D(trans->x, NULL, 0, xOut, NULL, xIn, yIn); + trans->y = psVectorFitPolynomial2D(trans->y, NULL, 0, yOut, NULL, xIn, yIn); + psFree(xIn); + psFree(yIn); + psFree(xOut); + psFree(yOut); + + if ((trans->x == NULL) || (trans->y == NULL)) { + psError( PS_ERR_UNKNOWN, true, "psVectorFitPolynomial2D() returned NULL: could not fit a 2-D polynomial to the data.\n"); + return(false); + } + + return(true); +} +