Index: /trunk/psLib/configure.ac =================================================================== --- /trunk/psLib/configure.ac (revision 3597) +++ /trunk/psLib/configure.ac (revision 3598) @@ -1,5 +1,5 @@ -AC_INIT([pslib],[1.4],[www.pan-starrs.org/bugzilla]) +AC_INIT([pslib],[1.5],[www.pan-starrs.org/bugzilla]) AM_CONFIG_HEADER(src/config.h) -AM_INIT_AUTOMAKE(pslib, 1.4) +AM_INIT_AUTOMAKE(pslib, 1.5) PSLIB_LT_VERSION="0:0:0" @@ -89,4 +89,5 @@ CFLAGS="${CFLAGS=} -DBUILD_PSDB" + PSLIB_CFLAGS="${PSLIB_CFLAGS=} -DBUILD_PSDB" else AC_MSG_RESULT([disable building MySQL functionality]) @@ -192,4 +193,6 @@ [PERL_PREFIX="$prefix"]) AC_SUBST(PERL_PREFIX,$PERL_PREFIX) + AC_SUBST(PERL_INSTALLSYTLE,[`perl '-V:installstyle'`]) + else AC_MSG_RESULT([disable building perl module]) Index: /trunk/psLib/pslib-config.in =================================================================== --- /trunk/psLib/pslib-config.in (revision 3597) +++ /trunk/psLib/pslib-config.in (revision 3598) @@ -4,4 +4,5 @@ exec_prefix=@exec_prefix@ includedir=@includedir@ +@PERL_INSTALLSYTLE@ usage() @@ -12,5 +13,6 @@ Known values for OPTION are: - --prefix show psLib installation prefix + --prefix print psLib installation prefix + --perlmodule print psLib perl module's installation location --libs print library linking information --cflags print pre-processor and compiler flags @@ -65,4 +67,7 @@ echo @PSLIB_LIBS@ ;; + --perlmodule) + echo @PERL_PREFIX@/$installstyle + ;; *) usage Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 3597) +++ /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); +} Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 3597) +++ /trunk/psLib/src/astro/psCoord.h (revision 3598) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-29 19:41:56 $ +* @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -344,8 +344,91 @@ ); -psSphere *psSpherePrecess(psSphere *coords, - const psTime *fromTime, - const psTime *toTime); - +/** Generates the complete spherical rotation to account for precession + * between two times. The equinoxes shall be Julian equinoxes. + * + * @return psSphere* the resulting spherical rotation + */ +psSphere* psSpherePrecess( + psSphere *coords, ///< coordinates (modified in-place) + const psTime *fromTime, ///< equinox of coords input + const psTime *toTime ///< equinox of coords output +); + +// XXX: Doxygenate. +psPlaneTransform *p_psPlaneTransformLinearInvert( + psPlaneTransform *transform +); + +// XXX: Doxygenate +psS32 p_psIsProjectionLinear( + psPlaneTransform *transform +); + +/** inverts a given transformation. + * + * It may assume that the input transformation is one-to-one, and that the + * inverse transformation may be specified through using polynomials of the + * same type and order as the forward transformation. In the event that the + * input transformation is linear, an exact solution may be calculated; + * otherwise nSamples samples in each axis, covering the region specified by + * region shall be used as a grid to fit the best inverse transformation. The + * function shall return NULL if it was unable to generate the inverse + * transformation; otherwise it shall return the inverse transformation. In + * the event that out is NULL, a new psPlaneTransform shall be allocated and + * returned. + * + * @return psPlaneTransform* the resulting inverted transform + */ +psPlaneTransform* psPlaneTransformInvert( + psPlaneTransform *out, ///< a transform to recycle, or NULL if one is to be created. + const psPlaneTransform *in, ///< transform to invert + psRegion *region, ///< region to fit for non-linear transform inversion + int nSamples ///< number of samples in each axis for fit +); + +/** Creates a single transformation that has the effect of performing trans1 + * followed by trans2. + * + * psPlaneTransformCombine takes two transformations (trans1 and trans2) and + * returns a single transformation that has the effect of performing trans1 + * followed by trans2. In the event that the input transformation is linear, + * an exact solution may be calculated; otherwise nSamples samples in each + * axis, covering the region specified by region shall be used as a grid to + * fit the best inverse transformation. The function shall return NULL if it + * was unable to generate the transformation; otherwise it shall return the + * transformation. + * + * @return psPlaneTransform* resulting transformation + */ +psPlaneTransform* psPlaneTransformCombine( + psPlaneTransform *out, ///< a transform to recycle, or NULL if one is to be created. + const psPlaneTransform *trans1, ///< first transform to combine + const psPlaneTransform *trans2 ///< first transform to combine +); + + +/** takes two arrays containing matched coordinates and returns the + * best-fitting transformation. + * + * psPlaneTransformFit takes two arrays containing matched coordinates (i.e., + * coordinates in the source array correspond to the coordinates in the dest + * array) and returns the best-fitting transformation. The source and dest + * will contain psCoords. In the event that the number of coordinates in each + * is not identical, the function shall generate a warning, and extra + * coordinates in the longer of the two shall be ignored. The trans transform + * may not be NULL, since it specifies the desired order, polynomial type and + * any polynomial terms to mask. nRejIter rejection iterations shall be + * performed, wherein coordinates lying more than sigmaClip standard + * deviations from the fit shall be rejected. + * + * @return bool TRUE if successful, otherwise FALSE. + */ +bool psPlaneTransformFit( + psPlaneTransform *trans, + const psArray *source, + const psArray *dest, + int nRejIter, + float sigmaClip +); /// @} Index: /trunk/psLib/src/astronomy/psAstrometry.c =================================================================== --- /trunk/psLib/src/astronomy/psAstrometry.c (revision 3597) +++ /trunk/psLib/src/astronomy/psAstrometry.c (revision 3598) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-30 02:21:14 $ + * @version $Revision: 1.63 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -51,106 +51,4 @@ } -/***************************************************************************** -isProjectionLinear(): 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. - *****************************************************************************/ -static psS32 isProjectionLinear(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); -} - -/***************************************************************************** -invertPlaneTransform(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. - - *****************************************************************************/ -static psPlaneTransform *invertPlaneTransform(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); -} static void FPAFree(psFPA* fpa) @@ -914,5 +812,5 @@ // generate an error if cell->toTP is not linear. - if (0 == isProjectionLinear(cell->toTP)) { + if (0 == p_psIsProjectionLinear(cell->toTP)) { psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psAstrometry_NONLINEAR_TRANSFORM, @@ -920,5 +818,5 @@ } - TPtoCell = invertPlaneTransform(cell->toTP); + TPtoCell = p_psPlaneTransformLinearInvert(cell->toTP); cellCoord = psPlaneTransformApply(cellCoord, TPtoCell, tpCoord); @@ -930,397 +828,3 @@ -/***************************************************************************** -multiplyCoeffs(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); -} - - - - -/***************************************************************************** -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 (isProjectionLinear((psPlaneTransform *) in)) { - printf("COOL: is linear\n"); - return(invertPlaneTransform((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); -} + Index: /trunk/psLib/src/astronomy/psAstrometry.h =================================================================== --- /trunk/psLib/src/astronomy/psAstrometry.h (revision 3597) +++ /trunk/psLib/src/astronomy/psAstrometry.h (revision 3598) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-31 02:33:28 $ +* @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -535,28 +535,3 @@ -// XXX: These functions don't belong here. Will migrate to psCoords.c later. -// XXX: Doxygenate. -psPlaneTransform *psPlaneTransformInvert( - psPlaneTransform *out, - const psPlaneTransform *in, - psRegion *region, - int nSamples -); - -// XXX: Doxygenate. -psPlaneTransform *psPlaneTransformCombine( - psPlaneTransform *out, - const psPlaneTransform *trans1, - const psPlaneTransform *trans2 -); - -// XXX: Doxygenate. -bool psPlaneTransformFit( - psPlaneTransform *trans, - const psArray *source, - const psArray *dest, - int nRejIter, - float sigmaClip -); - #endif Index: /trunk/psLib/src/astronomy/psCoord.c =================================================================== --- /trunk/psLib/src/astronomy/psCoord.c (revision 3597) +++ /trunk/psLib/src/astronomy/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); +} Index: /trunk/psLib/src/astronomy/psCoord.h =================================================================== --- /trunk/psLib/src/astronomy/psCoord.h (revision 3597) +++ /trunk/psLib/src/astronomy/psCoord.h (revision 3598) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-29 19:41:56 $ +* @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -344,8 +344,91 @@ ); -psSphere *psSpherePrecess(psSphere *coords, - const psTime *fromTime, - const psTime *toTime); - +/** Generates the complete spherical rotation to account for precession + * between two times. The equinoxes shall be Julian equinoxes. + * + * @return psSphere* the resulting spherical rotation + */ +psSphere* psSpherePrecess( + psSphere *coords, ///< coordinates (modified in-place) + const psTime *fromTime, ///< equinox of coords input + const psTime *toTime ///< equinox of coords output +); + +// XXX: Doxygenate. +psPlaneTransform *p_psPlaneTransformLinearInvert( + psPlaneTransform *transform +); + +// XXX: Doxygenate +psS32 p_psIsProjectionLinear( + psPlaneTransform *transform +); + +/** inverts a given transformation. + * + * It may assume that the input transformation is one-to-one, and that the + * inverse transformation may be specified through using polynomials of the + * same type and order as the forward transformation. In the event that the + * input transformation is linear, an exact solution may be calculated; + * otherwise nSamples samples in each axis, covering the region specified by + * region shall be used as a grid to fit the best inverse transformation. The + * function shall return NULL if it was unable to generate the inverse + * transformation; otherwise it shall return the inverse transformation. In + * the event that out is NULL, a new psPlaneTransform shall be allocated and + * returned. + * + * @return psPlaneTransform* the resulting inverted transform + */ +psPlaneTransform* psPlaneTransformInvert( + psPlaneTransform *out, ///< a transform to recycle, or NULL if one is to be created. + const psPlaneTransform *in, ///< transform to invert + psRegion *region, ///< region to fit for non-linear transform inversion + int nSamples ///< number of samples in each axis for fit +); + +/** Creates a single transformation that has the effect of performing trans1 + * followed by trans2. + * + * psPlaneTransformCombine takes two transformations (trans1 and trans2) and + * returns a single transformation that has the effect of performing trans1 + * followed by trans2. In the event that the input transformation is linear, + * an exact solution may be calculated; otherwise nSamples samples in each + * axis, covering the region specified by region shall be used as a grid to + * fit the best inverse transformation. The function shall return NULL if it + * was unable to generate the transformation; otherwise it shall return the + * transformation. + * + * @return psPlaneTransform* resulting transformation + */ +psPlaneTransform* psPlaneTransformCombine( + psPlaneTransform *out, ///< a transform to recycle, or NULL if one is to be created. + const psPlaneTransform *trans1, ///< first transform to combine + const psPlaneTransform *trans2 ///< first transform to combine +); + + +/** takes two arrays containing matched coordinates and returns the + * best-fitting transformation. + * + * psPlaneTransformFit takes two arrays containing matched coordinates (i.e., + * coordinates in the source array correspond to the coordinates in the dest + * array) and returns the best-fitting transformation. The source and dest + * will contain psCoords. In the event that the number of coordinates in each + * is not identical, the function shall generate a warning, and extra + * coordinates in the longer of the two shall be ignored. The trans transform + * may not be NULL, since it specifies the desired order, polynomial type and + * any polynomial terms to mask. nRejIter rejection iterations shall be + * performed, wherein coordinates lying more than sigmaClip standard + * deviations from the fit shall be rejected. + * + * @return bool TRUE if successful, otherwise FALSE. + */ +bool psPlaneTransformFit( + psPlaneTransform *trans, + const psArray *source, + const psArray *dest, + int nRejIter, + float sigmaClip +); /// @} Index: /trunk/psLib/src/astronomy/psDB.h =================================================================== --- /trunk/psLib/src/astronomy/psDB.h (revision 3597) +++ /trunk/psLib/src/astronomy/psDB.h (revision 3598) @@ -10,6 +10,6 @@ * @author Joshua Hoblitt * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-31 02:54:28 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2005 Joshua Hoblitt, University of Hawaii @@ -21,6 +21,4 @@ #ifdef BUILD_PSDB -#include - #include "psType.h" #include "psMetadata.h" @@ -36,5 +34,5 @@ typedef struct { - MYSQL *mysql; ///< MySQL database handle + void* mysql; ///< MySQL database handle } psDB; Index: /trunk/psLib/src/dataIO/psDB.h =================================================================== --- /trunk/psLib/src/dataIO/psDB.h (revision 3597) +++ /trunk/psLib/src/dataIO/psDB.h (revision 3598) @@ -10,6 +10,6 @@ * @author Joshua Hoblitt * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-31 02:54:28 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2005 Joshua Hoblitt, University of Hawaii @@ -21,6 +21,4 @@ #ifdef BUILD_PSDB -#include - #include "psType.h" #include "psMetadata.h" @@ -36,5 +34,5 @@ typedef struct { - MYSQL *mysql; ///< MySQL database handle + void* mysql; ///< MySQL database handle } psDB; Index: /trunk/psLib/src/dataManip/psConstants.h =================================================================== --- /trunk/psLib/src/dataManip/psConstants.h (revision 3597) +++ /trunk/psLib/src/dataManip/psConstants.h (revision 3598) @@ -1,3 +1,3 @@ -/** @file psComments.h +/** @file psConstants.h * * This file will hold definitions of various constants as well as common @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-30 23:44:05 $ + * @version $Revision: 1.63 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -13,6 +13,4 @@ * XXX: Add parenthesis around all arguments so that these macros can be * called with complex expressions. - * - * * */ Index: /trunk/psLib/src/dataManip/psFunctions.h =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.h (revision 3597) +++ /trunk/psLib/src/dataManip/psFunctions.h (revision 3598) @@ -1,20 +1,20 @@ /** @file psFunctions.h -* \brief Standard Mathematical Functions. -* \ingroup Stats -* -* This file will hold the prototypes for procedures which allocate, free, -* and evaluate various polynomials. Those polynomial structures are also -* defined here. -* -* @ingroup Stats -* -* @author Someone at IfA -* @author GLG, MHPCC -* -* @version $Revision: 1.43 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:23 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ + * @brief Standard Mathematical Functions. + * @ingroup Stats + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate various polynomials. Those polynomial structures are also + * defined here. + * + * @ingroup Stats + * + * @author Someone at IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ #if !defined(PS_FUNCTIONS_H) Index: /trunk/psLib/src/db/psDB.h =================================================================== --- /trunk/psLib/src/db/psDB.h (revision 3597) +++ /trunk/psLib/src/db/psDB.h (revision 3598) @@ -10,6 +10,6 @@ * @author Joshua Hoblitt * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-31 02:54:28 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2005 Joshua Hoblitt, University of Hawaii @@ -21,6 +21,4 @@ #ifdef BUILD_PSDB -#include - #include "psType.h" #include "psMetadata.h" @@ -36,5 +34,5 @@ typedef struct { - MYSQL *mysql; ///< MySQL database handle + void* mysql; ///< MySQL database handle } psDB; Index: /trunk/psLib/src/fileUtils/psDB.h =================================================================== --- /trunk/psLib/src/fileUtils/psDB.h (revision 3597) +++ /trunk/psLib/src/fileUtils/psDB.h (revision 3598) @@ -10,6 +10,6 @@ * @author Joshua Hoblitt * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-31 02:54:28 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2005 Joshua Hoblitt, University of Hawaii @@ -21,6 +21,4 @@ #ifdef BUILD_PSDB -#include - #include "psType.h" #include "psMetadata.h" @@ -36,5 +34,5 @@ typedef struct { - MYSQL *mysql; ///< MySQL database handle + void* mysql; ///< MySQL database handle } psDB; Index: /trunk/psLib/src/mainpage.dox =================================================================== --- /trunk/psLib/src/mainpage.dox (revision 3597) +++ /trunk/psLib/src/mainpage.dox (revision 3598) @@ -1,4 +1,14 @@ +/** @file mainpage.dox + * @brief Main page for the Doxygen documentation + * + * @author Robert Daniel DeSonia, MHPCC + * + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + /** @mainpage psLib Image Processing Library - @section intro Introduction Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 3597) +++ /trunk/psLib/src/math/psConstants.h (revision 3598) @@ -1,3 +1,3 @@ -/** @file psComments.h +/** @file psConstants.h * * This file will hold definitions of various constants as well as common @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-30 23:44:05 $ + * @version $Revision: 1.63 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -13,6 +13,4 @@ * XXX: Add parenthesis around all arguments so that these macros can be * called with complex expressions. - * - * * */ Index: /trunk/psLib/src/math/psPolynomial.h =================================================================== --- /trunk/psLib/src/math/psPolynomial.h (revision 3597) +++ /trunk/psLib/src/math/psPolynomial.h (revision 3598) @@ -1,20 +1,20 @@ /** @file psFunctions.h -* \brief Standard Mathematical Functions. -* \ingroup Stats -* -* This file will hold the prototypes for procedures which allocate, free, -* and evaluate various polynomials. Those polynomial structures are also -* defined here. -* -* @ingroup Stats -* -* @author Someone at IfA -* @author GLG, MHPCC -* -* @version $Revision: 1.43 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:23 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ + * @brief Standard Mathematical Functions. + * @ingroup Stats + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate various polynomials. Those polynomial structures are also + * defined here. + * + * @ingroup Stats + * + * @author Someone at IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ #if !defined(PS_FUNCTIONS_H) Index: /trunk/psLib/src/math/psSpline.h =================================================================== --- /trunk/psLib/src/math/psSpline.h (revision 3597) +++ /trunk/psLib/src/math/psSpline.h (revision 3598) @@ -1,20 +1,20 @@ /** @file psFunctions.h -* \brief Standard Mathematical Functions. -* \ingroup Stats -* -* This file will hold the prototypes for procedures which allocate, free, -* and evaluate various polynomials. Those polynomial structures are also -* defined here. -* -* @ingroup Stats -* -* @author Someone at IfA -* @author GLG, MHPCC -* -* @version $Revision: 1.43 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:23 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ + * @brief Standard Mathematical Functions. + * @ingroup Stats + * + * This file will hold the prototypes for procedures which allocate, free, + * and evaluate various polynomials. Those polynomial structures are also + * defined here. + * + * @ingroup Stats + * + * @author Someone at IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ #if !defined(PS_FUNCTIONS_H) Index: /trunk/psLib/src/psTest.h =================================================================== --- /trunk/psLib/src/psTest.h (revision 3597) +++ /trunk/psLib/src/psTest.h (revision 3598) @@ -1,2 +1,15 @@ +/** @file psTest.h + * + * Testing infrastructure functions. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + * + */ + #ifndef PSTEST_H Index: /trunk/psLib/src/pslib.h =================================================================== --- /trunk/psLib/src/pslib.h (revision 3597) +++ /trunk/psLib/src/pslib.h (revision 3598) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-31 03:02:15 $ +* @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -23,5 +23,4 @@ /******************************************************************************/ -#include "psDB.h" // // System Utilities @@ -169,4 +168,5 @@ #include "psAstrometry.h" #include "psConstants.h" +#include "psDB.h" /// @} Index: /trunk/psLib/swig/Makefile.am =================================================================== --- /trunk/psLib/swig/Makefile.am (revision 3597) +++ /trunk/psLib/swig/Makefile.am (revision 3598) @@ -23,5 +23,5 @@ psLibModule/Makefile.PL: psLibModule - cp $(srcdir)/Makefile.PL psLibModule/Makefile.PL + cp -f $(srcdir)/Makefile.PL psLibModule/Makefile.PL psLibModule/setup.txt:psLibModule Makefile @@ -50,3 +50,3 @@ uninstall-local: psLibModule/Makefile - cd PsLibModule && make clean + cd psLibModule && make clean Index: /trunk/psLib/test/astronomy/tst_psDB.c =================================================================== --- /trunk/psLib/test/astronomy/tst_psDB.c (revision 3597) +++ /trunk/psLib/test/astronomy/tst_psDB.c (revision 3598) @@ -9,12 +9,10 @@ * @author Aaron Culliney, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-30 23:22:39 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii * */ - -#define PS_ALLOW_MALLOC // for MySQL internals ... #include @@ -26,5 +24,4 @@ #include -#include "psDB.h" #include "psTest.h" #include "pslib.h" Index: /trunk/psLib/test/dataIO/tst_psDB.c =================================================================== --- /trunk/psLib/test/dataIO/tst_psDB.c (revision 3597) +++ /trunk/psLib/test/dataIO/tst_psDB.c (revision 3598) @@ -9,12 +9,10 @@ * @author Aaron Culliney, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-30 23:22:39 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii * */ - -#define PS_ALLOW_MALLOC // for MySQL internals ... #include @@ -26,5 +24,4 @@ #include -#include "psDB.h" #include "psTest.h" #include "pslib.h" Index: /trunk/psLib/test/dataIO/tst_psLookupTable_01.c =================================================================== --- /trunk/psLib/test/dataIO/tst_psLookupTable_01.c (revision 3597) +++ /trunk/psLib/test/dataIO/tst_psLookupTable_01.c (revision 3598) @@ -12,6 +12,6 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-25 19:55:01 $ +* @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii @@ -168,10 +168,10 @@ // Allocate lookup table with valid parameters - table1 = psLookupTableAlloc("verified/tableF32.dat", tableF32_validFrom, tableF32_validTo); + table1 = psLookupTableAlloc("tableF32.dat", tableF32_validFrom, tableF32_validTo); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Null lookup table generated from valid parameters"); return 1; } - if(strcmp(table1->fileName,"verified/tableF32.dat") != 0) { + if(strcmp(table1->fileName,"tableF32.dat") != 0) { psError(PS_ERR_UNKNOWN,true,"File name not properly stored in psLookupTable structure."); return 2; @@ -201,5 +201,5 @@ // Allocate table using table with psU8 index and valid types and index-values - table1 = psLookupTableAlloc("verified/tableU8.dat", 0, 100.5); + table1 = psLookupTableAlloc("tableU8.dat", 0, 100.5); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableU8.dat"); @@ -274,5 +274,5 @@ // Allocate table using table with psS32 index and valid types and index-values - table1 = psLookupTableAlloc("verified/tableS32.dat", -110, 1000.5); + table1 = psLookupTableAlloc("tableS32.dat", -110, 1000.5); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableS32.dat"); @@ -299,5 +299,5 @@ // Allocate table using table with psF32 index and valid types and index-values - table1 = psLookupTableAlloc("verified/tableF32.dat", -1100, 5500.5); + table1 = psLookupTableAlloc("tableF32.dat", -1100, 5500.5); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableF32.dat"); @@ -325,5 +325,5 @@ // Allocate table using psU8 index but with unsorted rows and verify the list is // sorted properly after being read - table1 = psLookupTableAlloc("verified/table10.dat",1,2); + table1 = psLookupTableAlloc("table10.dat",1,2); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table10.dat"); @@ -397,5 +397,5 @@ // Allocate table using table with invalid type in type row but valid index-values - table1 = psLookupTableAlloc("verified/table2.dat",0,99.99); + table1 = psLookupTableAlloc("table2.dat",0,99.99); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table2.dat"); @@ -411,5 +411,5 @@ // Allocate table using table with invalid value - table1 = psLookupTableAlloc("verified/table3.dat",0,75.0); + table1 = psLookupTableAlloc("table3.dat",0,75.0); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table3.dat"); @@ -464,5 +464,5 @@ // Interpolate values within the list and verify return values - table1 = psLookupTableAlloc("verified/table10.dat",0,20); + table1 = psLookupTableAlloc("table10.dat",0,20); table1 = psLookupTableRead(table1); for(psS32 i = 0; i < table1->numRows-1; i++ ) { @@ -547,5 +547,5 @@ // Interpolate values within the list and verify return values - table1 = psLookupTableAlloc("verified/table10.dat",0,20); + table1 = psLookupTableAlloc("table10.dat",0,20); table1 = psLookupTableRead(table1); statusVector = psVectorAlloc(table1->numCols,PS_TYPE_U32); Index: /trunk/psLib/test/fileUtils/tst_psDB.c =================================================================== --- /trunk/psLib/test/fileUtils/tst_psDB.c (revision 3597) +++ /trunk/psLib/test/fileUtils/tst_psDB.c (revision 3598) @@ -9,12 +9,10 @@ * @author Aaron Culliney, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-30 23:22:39 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii * */ - -#define PS_ALLOW_MALLOC // for MySQL internals ... #include @@ -26,5 +24,4 @@ #include -#include "psDB.h" #include "psTest.h" #include "pslib.h" Index: /trunk/psLib/test/fileUtils/tst_psLookupTable_01.c =================================================================== --- /trunk/psLib/test/fileUtils/tst_psLookupTable_01.c (revision 3597) +++ /trunk/psLib/test/fileUtils/tst_psLookupTable_01.c (revision 3598) @@ -12,6 +12,6 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-25 19:55:01 $ +* @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-31 23:01:46 $ * * Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii @@ -168,10 +168,10 @@ // Allocate lookup table with valid parameters - table1 = psLookupTableAlloc("verified/tableF32.dat", tableF32_validFrom, tableF32_validTo); + table1 = psLookupTableAlloc("tableF32.dat", tableF32_validFrom, tableF32_validTo); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Null lookup table generated from valid parameters"); return 1; } - if(strcmp(table1->fileName,"verified/tableF32.dat") != 0) { + if(strcmp(table1->fileName,"tableF32.dat") != 0) { psError(PS_ERR_UNKNOWN,true,"File name not properly stored in psLookupTable structure."); return 2; @@ -201,5 +201,5 @@ // Allocate table using table with psU8 index and valid types and index-values - table1 = psLookupTableAlloc("verified/tableU8.dat", 0, 100.5); + table1 = psLookupTableAlloc("tableU8.dat", 0, 100.5); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableU8.dat"); @@ -274,5 +274,5 @@ // Allocate table using table with psS32 index and valid types and index-values - table1 = psLookupTableAlloc("verified/tableS32.dat", -110, 1000.5); + table1 = psLookupTableAlloc("tableS32.dat", -110, 1000.5); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableS32.dat"); @@ -299,5 +299,5 @@ // Allocate table using table with psF32 index and valid types and index-values - table1 = psLookupTableAlloc("verified/tableF32.dat", -1100, 5500.5); + table1 = psLookupTableAlloc("tableF32.dat", -1100, 5500.5); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableF32.dat"); @@ -325,5 +325,5 @@ // Allocate table using psU8 index but with unsorted rows and verify the list is // sorted properly after being read - table1 = psLookupTableAlloc("verified/table10.dat",1,2); + table1 = psLookupTableAlloc("table10.dat",1,2); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table10.dat"); @@ -397,5 +397,5 @@ // Allocate table using table with invalid type in type row but valid index-values - table1 = psLookupTableAlloc("verified/table2.dat",0,99.99); + table1 = psLookupTableAlloc("table2.dat",0,99.99); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table2.dat"); @@ -411,5 +411,5 @@ // Allocate table using table with invalid value - table1 = psLookupTableAlloc("verified/table3.dat",0,75.0); + table1 = psLookupTableAlloc("table3.dat",0,75.0); if(table1 == NULL) { psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table3.dat"); @@ -464,5 +464,5 @@ // Interpolate values within the list and verify return values - table1 = psLookupTableAlloc("verified/table10.dat",0,20); + table1 = psLookupTableAlloc("table10.dat",0,20); table1 = psLookupTableRead(table1); for(psS32 i = 0; i < table1->numRows-1; i++ ) { @@ -547,5 +547,5 @@ // Interpolate values within the list and verify return values - table1 = psLookupTableAlloc("verified/table10.dat",0,20); + table1 = psLookupTableAlloc("table10.dat",0,20); table1 = psLookupTableRead(table1); statusVector = psVectorAlloc(table1->numCols,PS_TYPE_U32); Index: /trunk/psLib/test/sysUtils/verified/tst_psConfigure.stderr =================================================================== --- /trunk/psLib/test/sysUtils/verified/tst_psConfigure.stderr (revision 3597) +++ /trunk/psLib/test/sysUtils/verified/tst_psConfigure.stderr (revision 3598) @@ -6,5 +6,5 @@