Index: trunk/psLib/src/astro/psCoord.c =================================================================== --- trunk/psLib/src/astro/psCoord.c (revision 3540) +++ trunk/psLib/src/astro/psCoord.c (revision 3598) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.60 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-29 19:41:56 $ +* @version $Revision: 1.61 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -26,4 +26,6 @@ #include "psLogMsg.h" #include "psAstronomyErrors.h" +#include "psAstrometry.h" +#include "psMatrix.h" #include #include @@ -59,4 +61,39 @@ /* FUNCTION IMPLEMENTATION - LOCAL */ /*****************************************************************************/ +/***************************************************************************** +multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and +multiplies them. Basically, for each non-zero coeff in the trans1 coeff[][] +array, you must multiply by all non-zero coeffs in trans2. + +XXX: Inefficient in that the out polynomial is allocated every time. + *****************************************************************************/ + +psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1, + psDPolynomial2D *trans2) +{ + psS32 orderX = (trans1->nX + trans2->nX) - 1; + psS32 orderY = (trans1->nX + trans2->nX) - 1; + + psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD); + for (psS32 i = 0 ; i < out->nX; i++) { + for (psS32 j = 0 ; j < out->nY; j++) { + out->coeff[i][j] = 0.0; + out->mask[i][j] = 0; + } + } + + for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) { + for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) { + if (0.0 != trans1->coeff[t1x][t1y]) { + for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) { + for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) { + out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]); + } + } + } + } + } + return(out); +} /*****************************************************************************/ @@ -66,4 +103,107 @@ { // There are non dynamic allocated items +} + +/***************************************************************************** +p_psPlaneTransformLinearInvert(transform): : this is a private function which +simply inverts the supplied psPlaneTransform transform. It assumes that +"transform" is linear. + +This program assumes that the inverse of the following linear equations: + X2 = A + (B * X1) + (C * Y1); + Y2 = D + (E * X1) + (F * Y1); +is + Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B)); + X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C)); +or + X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) + + (X2 * -((F/C) / (E - ((F*B)/C)))) + + (Y2 * (1.0 / (E - ((F*B)/C)))); + Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) + + (X2 * -((E/B) / (F - ((C*E)/B)))) + + (Y2 * (1.0 / (F - ((C*E)/B)))); + +XXX: Since thre is now a general psPlaneTransformInvert() function, we +should rename this. + + *****************************************************************************/ +psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform) +{ + PS_PTR_CHECK_NULL(transform, 0); + PS_PTR_CHECK_NULL(transform->x, 0); + PS_PTR_CHECK_NULL(transform->y, 0); + + psF64 A = 0.0; + psF64 B = 0.0; + psF64 C = 0.0; + psF64 D = 0.0; + psF64 E = 0.0; + psF64 F = 0.0; + + // XXX: Test this for correctness. + A = transform->x->coeff[0][0]; + if (transform->x->nX >= 2) { + B = transform->x->coeff[1][0]; + } + if (transform->x->nY >= 2) { + C = transform->x->coeff[0][1]; + } + D = transform->y->coeff[0][0]; + if (transform->y->nX >= 2) { + E = transform->y->coeff[1][0]; + } + if (transform->y->nY >= 2) { + F = transform->y->coeff[0][1]; + } + + // XXX: Use the constructor here. + psPlaneTransform *out = psPlaneTransformAlloc(2, 2); + + out->x->coeff[0][0] = -D + ((F*A)/C) / (E - ((F*B)/C)); + out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C)); + out->x->coeff[0][1] = 1.0 / (E - ((F*B)/C)); + out->y->coeff[0][0] = -D + ((E*A)/B) / (F - ((C*E)/B)); + out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B)); + out->y->coeff[0][1] = 1.0 / (F - ((C*E)/B)); + + return(out); +} + +/***************************************************************************** +p_psIsProjectionLinear(): this is a private function which simply determines +if the supplied psPlaneTransform transform is linear: if any of the +cooefficients of order 2 are higher are non-zero, then it is not linear. + *****************************************************************************/ +psS32 p_psIsProjectionLinear(psPlaneTransform *transform) +{ + PS_PTR_CHECK_NULL(transform, 0); + PS_PTR_CHECK_NULL(transform->x, 0); + PS_PTR_CHECK_NULL(transform->y, 0); + + for (psS32 i=0;i<(transform->x->nX);i++) { + for (psS32 j=0;j<(transform->x->nY);j++) { + if (transform->x->coeff[i][j] != 0.0) { + if (!(((i == 0) && (j == 0)) || + ((i == 0) && (j == 1)) || + ((i == 1) && (j == 0)))) { + return(0); + } + } + } + } + + for (psS32 i=0;i<(transform->y->nX);i++) { + for (psS32 j=0;j<(transform->y->nY);j++) { + if (transform->y->coeff[i][j] != 0.0) { + if (!(((i == 0) && (j == 0)) || + ((i == 0) && (j == 1)) || + ((i == 1) && (j == 0)))) { + return(0); + } + } + } + } + + return(1); } @@ -721,2 +861,361 @@ } + + +/***************************************************************************** +psPlaneTransformCombine(out, trans1, trans2) + +XXX: Much room for optimization. Currently, we call the polyMultiply +routine far too many times. + *****************************************************************************/ +psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out, + const psPlaneTransform *trans1, + const psPlaneTransform *trans2) +{ + PS_PTR_CHECK_NULL(trans1, NULL); + PS_PTR_CHECK_NULL(trans2, NULL); + + // + // Determine the size of the new psPlaneTransform. + // + // PS_MAX( Number of x terms in T2->x * number of x terms in T1->x, + // Number of y terms in T2->x * number of x terms in T1->y, + psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX), + (trans2->x->nY * trans1->y->nX)); + psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY), + (trans2->x->nY * trans1->y->nY)); + + psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX), + (trans2->y->nY * trans1->y->nX)); + psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY), + (trans2->y->nY * trans1->y->nY)); + psS32 orderX = PS_MAX(orderXnX, orderYnX); + psS32 orderY = PS_MAX(orderXnY, orderYnY); + + // + // Allocate the new psPlaneTransform, if necessary. + // + psPlaneTransform *myPT = NULL; + if (out == NULL) { + myPT = psPlaneTransformAlloc(orderX, orderY); + } else { + if ((out->x->nX == orderX) && (out->x->nY == orderY) && + (out->y->nX == orderX) && (out->y->nY == orderY)) { + myPT = out; + } else { + psFree(out); + myPT = psPlaneTransformAlloc(orderX, orderY); + } + } + + // + // Initialize the new psPlaneTransform, if necessary. + // + for (psS32 i = 0 ; i < orderX ; i++) { + for (psS32 j = 0 ; j < orderY ; j++) { + myPT->x->coeff[i][j] = 0.0; + myPT->x->mask[i][j] = 0; + myPT->y->coeff[i][j] = 0.0; + myPT->y->mask[i][j] = 0; + } + } + + // + // For each term (a * x^i * y^j) in trans2, we substitute the appropriate + // equation from trans1, and raise it to the appropriate power. This is + // done via the multiplyDPoly2D(). The result is a polynomial (currPoly) + // and its coefficients are added into the myPT coeff matrix. + // + // XXX: This is horribly inefficient in that the trans1 polys are repeatedly + // multiplied against themselves. This can easily be improved. + // + for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) { + for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) { + psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); + currPoly->coeff[0][0] = 1.0; + currPoly->mask[0][0] = 0; + psDPolynomial2D *newPoly = NULL; + + if (trans2->x->mask[t2x][t2y] == 0) { + + // Must raise trans1->y to the t2y-power. + for (psS32 c = 0 ; c < t2y; c++) { + newPoly = multiplyDPoly2D(currPoly, trans1->y); + psFree(currPoly); + currPoly = newPoly; + } + + // Must raise trans1->x to the t2x-power. + for (psS32 c = 0 ; c < t2x; c++) { + newPoly = multiplyDPoly2D(currPoly, trans1->x); + psFree(currPoly); + currPoly = newPoly; + } + + // Set the appropriate coeffs in myPT->x + for (psS32 i = 0 ; i < currPoly->nX ; i++) { + for (psS32 j = 0 ; j < currPoly->nY ; j++) { + myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y]; + } + } + } + psFree(currPoly); + } + } + + + + for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) { + for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) { + psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); + currPoly->coeff[0][0] = 1.0; + currPoly->mask[0][0] = 0; + psDPolynomial2D *newPoly = NULL; + + if (trans2->y->mask[t2x][t2y] == 0) { + + // Must raise trans1->y to the t2y-power. + for (psS32 c = 0 ; c < t2y; c++) { + newPoly = multiplyDPoly2D(currPoly, trans1->y); + psFree(currPoly); + currPoly = newPoly; + } + + // Must raise trans1->x to the t2x-power. + for (psS32 c = 0 ; c < t2x; c++) { + newPoly = multiplyDPoly2D(currPoly, trans1->x); + psFree(currPoly); + currPoly = newPoly; + } + + // Set the appropriate coeffs in myPT->x + for (psS32 i = 0 ; i < currPoly->nX ; i++) { + for (psS32 j = 0 ; j < currPoly->nY ; j++) { + myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y]; + } + } + } + psFree(currPoly); + } + } + + return(myPT); +} + +/***************************************************************************** +psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip) + +XXX: What about nRejIter? Iterations? +XXX: Use static vectors for internal data. + *****************************************************************************/ +bool psPlaneTransformFit(psPlaneTransform *trans, + const psArray *source, + const psArray *dest, + int nRejIter, + float sigmaClip) +{ + PS_PTR_CHECK_NULL(trans, NULL); + PS_PTR_CHECK_NULL(source, NULL); + PS_PTR_CHECK_NULL(dest, NULL); + + psS32 numCoords = PS_MIN(source->n, dest->n); + // This is not really necessary because of above conditionals. + psS32 order = PS_MAX(trans->x->nX, trans->x->nY); + + // + // Create fake polynomial to use in evaluation + // + psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD); + 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 + // + 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 = psDPolynomial2DEval(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 = psDPolynomial2DEval(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 *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition + psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix + psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x + psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y + + // + // XXX: Should check the output of the matrix routines and return false if bad. + // + + // + // 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(true); +} + + +/***************************************************************************** +psPlaneTransformInvert(out, in, region, nSamples) + +// XXX: Use static data structures. + *****************************************************************************/ +psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out, + const psPlaneTransform *in, + psRegion *region, + int nSamples) +{ + PS_PTR_CHECK_NULL(in, NULL); + // + // If the transform is linear, then invert it exactly and return. + // + if (p_psIsProjectionLinear((psPlaneTransform *) in)) { + printf("COOL: is linear\n"); + return(p_psPlaneTransformLinearInvert((psPlaneTransform *) in)); + } + PS_PTR_CHECK_NULL(region, NULL); + PS_INT_COMPARE(1, nSamples, NULL); + + // Ensure that the input transformation is symmetrical. + if ((in->x->nX != in->x->nY) || + (in->y->nX != in->y->nY) || + (in->x->nX != in->y->nX)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY."); + } + psS32 order = PS_MAX(in->x->nX, in->x->nY); + + psPlaneTransform *myPT = NULL; + psPlane *inCoord = psPlaneAlloc(); + psPlane *outCoord = psPlaneAlloc(); + + // + // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size. + // + if (out == NULL) { + myPT = psPlaneTransformAlloc(order, order); + } else { + if ((out->x->nX == order) && (out->x->nY == order) && + (out->y->nX == order) && (out->y->nY == order)) { + myPT = out; + } else { + psFree(out); + myPT = psPlaneTransformAlloc(order, order); + } + } + + // + // Copy the input transform to myPT. + // + for (psS32 i = 0 ; i < in->x->nX ; i++) { + for (psS32 j = 0 ; j < in->x->nY ; j++) { + myPT->x->coeff[i][j] = in->x->coeff[i][j]; + } + } + for (psS32 i = 0 ; i < in->y->nX ; i++) { + for (psS32 j = 0 ; j < in->y->nY ; j++) { + myPT->y->coeff[i][j] = in->y->coeff[i][j]; + } + } + + // + // Create a grid of xin,yin --> xout,yout + // + psArray *inData = psArrayAlloc(nSamples * nSamples); + psArray *outData = psArrayAlloc(nSamples * nSamples); + for (psS32 i = 0 ; i < inData->n; i++) { + inData->data[i] = (psPtr *) psPlaneAlloc(); + outData->data[i] = (psPtr *) psPlaneAlloc(); + } + + // + // Initialize the grid. + // + psS32 cnt = 0; + for (int yint = 0; yint < nSamples; yint++) { + inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples)); + for (int xint = 0; xint < nSamples; xint++) { + inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples)); + (void)psPlaneTransformApply(outCoord, in, inCoord); + + ((psPlane *) outData->data[cnt])->x = inCoord->x; + ((psPlane *) outData->data[cnt])->y = inCoord->y; + ((psPlane *) inData->data[cnt])->x = outCoord->x; + ((psPlane *) inData->data[cnt])->y = outCoord->y; + + cnt++; + } + } + bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0); + + psFree(inCoord); + psFree(outCoord); + psFree(inData); + psFree(outData); + + if (rc == true) { + return(myPT); + } + + // XXX: Generate an error message, or warning message. + return(NULL); +}