IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3598


Ignore:
Timestamp:
Mar 31, 2005, 1:01:46 PM (21 years ago)
Author:
desonia
Message:

cosmetic tweaks.

Location:
trunk/psLib
Files:
27 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/configure.ac

    r3588 r3598  
    1 AC_INIT([pslib],[1.4],[www.pan-starrs.org/bugzilla])
     1AC_INIT([pslib],[1.5],[www.pan-starrs.org/bugzilla])
    22AM_CONFIG_HEADER(src/config.h)
    3 AM_INIT_AUTOMAKE(pslib, 1.4)
     3AM_INIT_AUTOMAKE(pslib, 1.5)
    44
    55PSLIB_LT_VERSION="0:0:0"
    … …  
    8989
    9090  CFLAGS="${CFLAGS=} -DBUILD_PSDB"
     91  PSLIB_CFLAGS="${PSLIB_CFLAGS=} -DBUILD_PSDB"
    9192else
    9293  AC_MSG_RESULT([disable building MySQL functionality])
    … …  
    192193    [PERL_PREFIX="$prefix"])
    193194    AC_SUBST(PERL_PREFIX,$PERL_PREFIX)
     195    AC_SUBST(PERL_INSTALLSYTLE,[`perl '-V:installstyle'`])
     196
    194197else
    195198  AC_MSG_RESULT([disable building perl module])

  • trunk/psLib/pslib-config.in

    r3237 r3598  
    44exec_prefix=@exec_prefix@
    55includedir=@includedir@
     6@PERL_INSTALLSYTLE@
    67
    78usage()
    … …  
    1213Known values for OPTION are:
    1314
    14   --prefix              show psLib installation prefix
     15  --prefix              print psLib installation prefix
     16  --perlmodule          print psLib perl module's installation location
    1517  --libs                print library linking information
    1618  --cflags              print pre-processor and compiler flags
    … …  
    6567        echo @PSLIB_LIBS@
    6668        ;;
     69    --perlmodule)
     70        echo @PERL_PREFIX@/$installstyle
     71        ;;
    6772    *)
    6873        usage

  • trunk/psLib/src/astro/psCoord.c

    r3540 r3598  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-29 19:41:56 $
     12*  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-31 23:01:46 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2626#include "psLogMsg.h"
    2727#include "psAstronomyErrors.h"
     28#include "psAstrometry.h"
     29#include "psMatrix.h"
    2830#include <math.h>
    2931#include <float.h>
    … …  
    5961/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
    6062/*****************************************************************************/
     63/*****************************************************************************
     64multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and
     65multiplies them.  Basically, for each non-zero coeff in the trans1 coeff[][]
     66array, you must multiply by all non-zero coeffs in trans2.
     67 
     68XXX: Inefficient in that the out polynomial is allocated every time.
     69 *****************************************************************************/
     70
     71psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1,
     72                                 psDPolynomial2D *trans2)
     73{
     74    psS32 orderX = (trans1->nX + trans2->nX) - 1;
     75    psS32 orderY = (trans1->nX + trans2->nX) - 1;
     76
     77    psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
     78    for (psS32 i = 0 ; i < out->nX; i++) {
     79        for (psS32 j = 0 ; j < out->nY; j++) {
     80            out->coeff[i][j] = 0.0;
     81            out->mask[i][j] = 0;
     82        }
     83    }
     84
     85    for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) {
     86        for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) {
     87            if (0.0 != trans1->coeff[t1x][t1y]) {
     88                for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
     89                    for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
     90                        out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
     91                    }
     92                }
     93            }
     94        }
     95    }
     96    return(out);
     97}
    6198
    6299/*****************************************************************************/
    … …  
    66103{
    67104    // There are non dynamic allocated items
     105}
     106
     107/*****************************************************************************
     108p_psPlaneTransformLinearInvert(transform): : this is a private function which
     109simply inverts the supplied psPlaneTransform transform.  It assumes that
     110"transform" is linear.
     111 
     112This program assumes that the inverse of the following linear equations:
     113        X2 = A + (B * X1) + (C * Y1);
     114        Y2 = D + (E * X1) + (F * Y1);
     115is
     116        Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B));
     117        X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
     118or
     119 X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) +
     120      (X2 * -((F/C) / (E - ((F*B)/C)))) +
     121             (Y2 * (1.0 / (E - ((F*B)/C))));
     122        Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) +
     123             (X2 * -((E/B) / (F - ((C*E)/B)))) +
     124             (Y2 * (1.0 / (F - ((C*E)/B))));
     125 
     126XXX: Since thre is now a general psPlaneTransformInvert() function, we
     127should rename this.
     128 
     129 *****************************************************************************/
     130psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
     131{
     132    PS_PTR_CHECK_NULL(transform, 0);
     133    PS_PTR_CHECK_NULL(transform->x, 0);
     134    PS_PTR_CHECK_NULL(transform->y, 0);
     135
     136    psF64 A = 0.0;
     137    psF64 B = 0.0;
     138    psF64 C = 0.0;
     139    psF64 D = 0.0;
     140    psF64 E = 0.0;
     141    psF64 F = 0.0;
     142
     143    // XXX: Test this for correctness.
     144    A = transform->x->coeff[0][0];
     145    if (transform->x->nX >= 2) {
     146        B = transform->x->coeff[1][0];
     147    }
     148    if (transform->x->nY >= 2) {
     149        C = transform->x->coeff[0][1];
     150    }
     151    D = transform->y->coeff[0][0];
     152    if (transform->y->nX >= 2) {
     153        E = transform->y->coeff[1][0];
     154    }
     155    if (transform->y->nY >= 2) {
     156        F = transform->y->coeff[0][1];
     157    }
     158
     159    // XXX: Use the constructor here.
     160    psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
     161
     162    out->x->coeff[0][0] = -D + ((F*A)/C) / (E - ((F*B)/C));
     163    out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C));
     164    out->x->coeff[0][1] =  1.0 / (E - ((F*B)/C));
     165    out->y->coeff[0][0] = -D + ((E*A)/B) / (F - ((C*E)/B));
     166    out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B));
     167    out->y->coeff[0][1] =  1.0 / (F - ((C*E)/B));
     168
     169    return(out);
     170}
     171
     172/*****************************************************************************
     173p_psIsProjectionLinear(): this is a private function which simply determines
     174if the supplied psPlaneTransform transform is linear: if any of the
     175cooefficients of order 2 are higher are non-zero, then it is not linear.
     176 *****************************************************************************/
     177psS32 p_psIsProjectionLinear(psPlaneTransform *transform)
     178{
     179    PS_PTR_CHECK_NULL(transform, 0);
     180    PS_PTR_CHECK_NULL(transform->x, 0);
     181    PS_PTR_CHECK_NULL(transform->y, 0);
     182
     183    for (psS32 i=0;i<(transform->x->nX);i++) {
     184        for (psS32 j=0;j<(transform->x->nY);j++) {
     185            if (transform->x->coeff[i][j] != 0.0) {
     186                if (!(((i == 0) && (j == 0)) ||
     187                        ((i == 0) && (j == 1)) ||
     188                        ((i == 1) && (j == 0)))) {
     189                    return(0);
     190                }
     191            }
     192        }
     193    }
     194
     195    for (psS32 i=0;i<(transform->y->nX);i++) {
     196        for (psS32 j=0;j<(transform->y->nY);j++) {
     197            if (transform->y->coeff[i][j] != 0.0) {
     198                if (!(((i == 0) && (j == 0)) ||
     199                        ((i == 0) && (j == 1)) ||
     200                        ((i == 1) && (j == 0)))) {
     201                    return(0);
     202                }
     203            }
     204        }
     205    }
     206
     207    return(1);
    68208}
    69209
    … …  
    721861}
    722862
     863
     864
     865/*****************************************************************************
     866psPlaneTransformCombine(out, trans1, trans2)
     867 
     868XXX: Much room for optimization.  Currently, we call the polyMultiply
     869routine far too many times.
     870 *****************************************************************************/
     871psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out,
     872        const psPlaneTransform *trans1,
     873        const psPlaneTransform *trans2)
     874{
     875    PS_PTR_CHECK_NULL(trans1, NULL);
     876    PS_PTR_CHECK_NULL(trans2, NULL);
     877
     878    //
     879    // Determine the size of the new psPlaneTransform.
     880    //
     881    // PS_MAX(  Number of x terms in T2->x * number of x terms in T1->x,
     882    //          Number of y terms in T2->x * number of x terms in T1->y,
     883    psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX),
     884                            (trans2->x->nY * trans1->y->nX));
     885    psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY),
     886                            (trans2->x->nY * trans1->y->nY));
     887
     888    psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX),
     889                            (trans2->y->nY * trans1->y->nX));
     890    psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY),
     891                            (trans2->y->nY * trans1->y->nY));
     892    psS32 orderX = PS_MAX(orderXnX, orderYnX);
     893    psS32 orderY = PS_MAX(orderXnY, orderYnY);
     894
     895    //
     896    // Allocate the new psPlaneTransform, if necessary.
     897    //
     898    psPlaneTransform *myPT = NULL;
     899    if (out == NULL) {
     900        myPT = psPlaneTransformAlloc(orderX, orderY);
     901    } else {
     902        if ((out->x->nX == orderX) && (out->x->nY == orderY) &&
     903                (out->y->nX == orderX) && (out->y->nY == orderY)) {
     904            myPT = out;
     905        } else {
     906            psFree(out);
     907            myPT = psPlaneTransformAlloc(orderX, orderY);
     908        }
     909    }
     910
     911    //
     912    // Initialize the new psPlaneTransform, if necessary.
     913    //
     914    for (psS32 i = 0 ; i < orderX ; i++) {
     915        for (psS32 j = 0 ; j < orderY ; j++) {
     916            myPT->x->coeff[i][j] = 0.0;
     917            myPT->x->mask[i][j] = 0;
     918            myPT->y->coeff[i][j] = 0.0;
     919            myPT->y->mask[i][j] = 0;
     920        }
     921    }
     922
     923    //
     924    // For each term (a * x^i * y^j) in trans2, we substitute the appropriate
     925    // equation from trans1, and raise it to the appropriate power.  This is
     926    // done via the multiplyDPoly2D().  The result is a polynomial (currPoly)
     927    // and its coefficients are added into the myPT coeff matrix.
     928    //
     929    // XXX: This is horribly inefficient in that the trans1 polys are repeatedly
     930    // multiplied against themselves.  This can easily be improved.
     931    //
     932    for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) {
     933        for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) {
     934            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
     935            currPoly->coeff[0][0] = 1.0;
     936            currPoly->mask[0][0] = 0;
     937            psDPolynomial2D *newPoly = NULL;
     938
     939            if (trans2->x->mask[t2x][t2y] == 0) {
     940
     941                // Must raise trans1->y to the t2y-power.
     942                for (psS32 c = 0 ; c < t2y; c++) {
     943                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
     944                    psFree(currPoly);
     945                    currPoly = newPoly;
     946                }
     947
     948                // Must raise trans1->x to the t2x-power.
     949                for (psS32 c = 0 ; c < t2x; c++) {
     950                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
     951                    psFree(currPoly);
     952                    currPoly = newPoly;
     953                }
     954
     955                // Set the appropriate coeffs in myPT->x
     956                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
     957                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
     958                        myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
     959                    }
     960                }
     961            }
     962            psFree(currPoly);
     963        }
     964    }
     965
     966
     967
     968    for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) {
     969        for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) {
     970            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
     971            currPoly->coeff[0][0] = 1.0;
     972            currPoly->mask[0][0] = 0;
     973            psDPolynomial2D *newPoly = NULL;
     974
     975            if (trans2->y->mask[t2x][t2y] == 0) {
     976
     977                // Must raise trans1->y to the t2y-power.
     978                for (psS32 c = 0 ; c < t2y; c++) {
     979                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
     980                    psFree(currPoly);
     981                    currPoly = newPoly;
     982                }
     983
     984                // Must raise trans1->x to the t2x-power.
     985                for (psS32 c = 0 ; c < t2x; c++) {
     986                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
     987                    psFree(currPoly);
     988                    currPoly = newPoly;
     989                }
     990
     991                // Set the appropriate coeffs in myPT->x
     992                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
     993                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
     994                        myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
     995                    }
     996                }
     997            }
     998            psFree(currPoly);
     999        }
     1000    }
     1001
     1002    return(myPT);
     1003}
     1004
     1005/*****************************************************************************
     1006psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip)
     1007 
     1008XXX: What about nRejIter?  Iterations?
     1009XXX: Use static vectors for internal data.
     1010 *****************************************************************************/
     1011bool psPlaneTransformFit(psPlaneTransform *trans,
     1012                         const psArray *source,
     1013                         const psArray *dest,
     1014                         int nRejIter,
     1015                         float sigmaClip)
     1016{
     1017    PS_PTR_CHECK_NULL(trans, NULL);
     1018    PS_PTR_CHECK_NULL(source, NULL);
     1019    PS_PTR_CHECK_NULL(dest, NULL);
     1020
     1021    psS32 numCoords = PS_MIN(source->n, dest->n);
     1022    // This is not really necessary because of above conditionals.
     1023    psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
     1024
     1025    //
     1026    // Create fake polynomial to use in evaluation
     1027    //
     1028    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
     1029    for (int i = 0; i < order; i++) {
     1030        for (int j = 0; j < order; j++) {
     1031            fakePoly->coeff[i][j] = 1.0;
     1032            fakePoly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
     1033        }
     1034    }
     1035
     1036    //
     1037    // Initialize the matrix and vectors
     1038    //
     1039    psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
     1040    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
     1041    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
     1042    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
     1043    for (psS32 i = 0; i < nCoeff; i++) {
     1044        for (psS32 j = 0; j < nCoeff; j++) {
     1045            matrix->data.F64[i][j] = 0.0;
     1046        }
     1047        xVector->data.F64[i] = 0.0;
     1048        yVector->data.F64[i] = 0.0;
     1049    }
     1050
     1051    //
     1052    // Iterate over the grid points
     1053    //
     1054    for (psS32 g = 0; g < numCoords; g++) {
     1055        // Iterate over the polynomial coefficients, accumulating the matrix and vectors
     1056
     1057        for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1058            for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1059                fakePoly->mask[i][j] = 0;
     1060                psF64 xIn = ((psPlane *) source->data[g])->x;
     1061                psF64 yIn = ((psPlane *) source->data[g])->y;
     1062                psF64 xOut = ((psPlane *) dest->data[g])->x;
     1063                psF64 yOut = ((psPlane *) dest->data[g])->y;
     1064                psF64 ijPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
     1065                fakePoly->mask[i][j] = 1;
     1066
     1067                for (psS32 m = 0, mnIndex = 0; m < order; m++) {
     1068                    for (psS32 n = 0; n < order - m; n++, mnIndex++) {
     1069                        fakePoly->mask[m][n] = 0;
     1070                        psF64 mnPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
     1071                        fakePoly->mask[m][n] = 1;
     1072
     1073                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
     1074                    }
     1075                }
     1076
     1077                xVector->data.F64[ijIndex] += ijPoly * xOut;
     1078                yVector->data.F64[ijIndex] += ijPoly * yOut;
     1079            }
     1080        }
     1081    }
     1082
     1083    //
     1084    // Solution via LU Decomposition
     1085    //
     1086    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     1087    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
     1088    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
     1089    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
     1090
     1091    //
     1092    // XXX: Should check the output of the matrix routines and return false if bad.
     1093    //
     1094
     1095    //
     1096    // Stuff coefficients into transformation
     1097    //
     1098    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1099        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1100            trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
     1101            trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
     1102        }
     1103    }
     1104
     1105    psFree(fakePoly);
     1106    psFree(permutation);
     1107    psFree(luMatrix);
     1108    psFree(xSolution);
     1109    psFree(ySolution);
     1110    psFree(matrix);
     1111    psFree(xVector);
     1112    psFree(yVector);
     1113
     1114    return(true);
     1115}
     1116
     1117
     1118/*****************************************************************************
     1119psPlaneTransformInvert(out, in, region, nSamples)
     1120 
     1121// XXX: Use static data structures.
     1122 *****************************************************************************/
     1123psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
     1124        const psPlaneTransform *in,
     1125        psRegion *region,
     1126        int nSamples)
     1127{
     1128    PS_PTR_CHECK_NULL(in, NULL);
     1129    //
     1130    // If the transform is linear, then invert it exactly and return.
     1131    //
     1132    if (p_psIsProjectionLinear((psPlaneTransform *) in)) {
     1133        printf("COOL: is linear\n");
     1134        return(p_psPlaneTransformLinearInvert((psPlaneTransform *) in));
     1135    }
     1136    PS_PTR_CHECK_NULL(region, NULL);
     1137    PS_INT_COMPARE(1, nSamples, NULL);
     1138
     1139    // Ensure that the input transformation is symmetrical.
     1140    if ((in->x->nX != in->x->nY) ||
     1141            (in->y->nX != in->y->nY) ||
     1142            (in->x->nX != in->y->nX)) {
     1143        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
     1144    }
     1145    psS32 order = PS_MAX(in->x->nX, in->x->nY);
     1146
     1147    psPlaneTransform *myPT = NULL;
     1148    psPlane *inCoord = psPlaneAlloc();
     1149    psPlane *outCoord = psPlaneAlloc();
     1150
     1151    //
     1152    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
     1153    //
     1154    if (out == NULL) {
     1155        myPT = psPlaneTransformAlloc(order, order);
     1156    } else {
     1157        if ((out->x->nX == order) && (out->x->nY == order) &&
     1158                (out->y->nX == order) && (out->y->nY == order)) {
     1159            myPT = out;
     1160        } else {
     1161            psFree(out);
     1162            myPT = psPlaneTransformAlloc(order, order);
     1163        }
     1164    }
     1165
     1166    //
     1167    // Copy the input transform to myPT.
     1168    //
     1169    for (psS32 i = 0 ; i < in->x->nX ; i++) {
     1170        for (psS32 j = 0 ; j < in->x->nY ; j++) {
     1171            myPT->x->coeff[i][j] = in->x->coeff[i][j];
     1172        }
     1173    }
     1174    for (psS32 i = 0 ; i < in->y->nX ; i++) {
     1175        for (psS32 j = 0 ; j < in->y->nY ; j++) {
     1176            myPT->y->coeff[i][j] = in->y->coeff[i][j];
     1177        }
     1178    }
     1179
     1180    //
     1181    // Create a grid of xin,yin --> xout,yout
     1182    //
     1183    psArray *inData = psArrayAlloc(nSamples * nSamples);
     1184    psArray *outData = psArrayAlloc(nSamples * nSamples);
     1185    for (psS32 i = 0 ; i < inData->n; i++) {
     1186        inData->data[i] = (psPtr *) psPlaneAlloc();
     1187        outData->data[i] = (psPtr *) psPlaneAlloc();
     1188    }
     1189
     1190    //
     1191    // Initialize the grid.
     1192    //
     1193    psS32 cnt = 0;
     1194    for (int yint = 0; yint < nSamples; yint++) {
     1195        inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
     1196        for (int xint = 0; xint < nSamples; xint++) {
     1197            inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
     1198            (void)psPlaneTransformApply(outCoord, in, inCoord);
     1199
     1200            ((psPlane *) outData->data[cnt])->x = inCoord->x;
     1201            ((psPlane *) outData->data[cnt])->y = inCoord->y;
     1202            ((psPlane *) inData->data[cnt])->x = outCoord->x;
     1203            ((psPlane *) inData->data[cnt])->y = outCoord->y;
     1204
     1205            cnt++;
     1206        }
     1207    }
     1208    bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0);
     1209
     1210    psFree(inCoord);
     1211    psFree(outCoord);
     1212    psFree(inData);
     1213    psFree(outData);
     1214
     1215    if (rc == true) {
     1216        return(myPT);
     1217    }
     1218
     1219    // XXX: Generate an error message, or warning message.
     1220    return(NULL);
     1221}

  • trunk/psLib/src/astro/psCoord.h

    r3540 r3598  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-29 19:41:56 $
     12*  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-31 23:01:46 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    344344);
    345345
    346 psSphere *psSpherePrecess(psSphere *coords,
    347                           const psTime *fromTime,
    348                           const psTime *toTime);
    349 
     346/** Generates the complete spherical rotation to account for precession
     347 *  between two times.  The equinoxes shall be Julian equinoxes.
     348 *
     349 *  @return psSphere* the resulting spherical rotation
     350 */
     351psSphere* psSpherePrecess(
     352    psSphere *coords,                  ///< coordinates (modified in-place)
     353    const psTime *fromTime,            ///< equinox of coords input
     354    const psTime *toTime               ///< equinox of coords output
     355);
     356
     357// XXX: Doxygenate.
     358psPlaneTransform *p_psPlaneTransformLinearInvert(
     359    psPlaneTransform *transform
     360);
     361
     362// XXX: Doxygenate
     363psS32 p_psIsProjectionLinear(
     364    psPlaneTransform *transform
     365);
     366
     367/** inverts a given transformation.
     368 *
     369 *  It may assume that the input transformation is one-to-one, and that the
     370 *  inverse transformation may be specified through using polynomials of the
     371 *  same type and order as the forward transformation. In the event that the
     372 *  input transformation is linear, an exact solution may be calculated;
     373 *  otherwise nSamples samples in each axis, covering the region specified by
     374 *  region shall be used as a grid to fit the best inverse transformation. The
     375 *  function shall return NULL if it was unable to generate the inverse
     376 *  transformation; otherwise it shall return the inverse transformation. In
     377 *  the event that out is NULL, a new psPlaneTransform shall be allocated and
     378 *  returned.
     379 *
     380 *  @return psPlaneTransform*  the resulting inverted transform
     381 */
     382psPlaneTransform* psPlaneTransformInvert(
     383    psPlaneTransform *out,             ///< a transform to recycle, or NULL if one is to be created.
     384    const psPlaneTransform *in,        ///< transform to invert
     385    psRegion *region,                  ///< region to fit for non-linear transform inversion
     386    int nSamples                       ///< number of samples in each axis for fit
     387);
     388
     389/** Creates a single transformation that has the effect of performing trans1
     390 *  followed by trans2.
     391 *
     392 *  psPlaneTransformCombine takes two transformations (trans1 and trans2) and
     393 *  returns a single transformation that has the effect of performing trans1
     394 *  followed by trans2. In the event that the input transformation is linear,
     395 *  an exact solution may be calculated; otherwise nSamples samples in each
     396 *  axis, covering the region specified by region shall be used as a grid to
     397 *  fit the best inverse transformation. The function shall return NULL if it
     398 *  was unable to generate the transformation; otherwise it shall return the
     399 *  transformation.
     400 *
     401 *  @return psPlaneTransform*    resulting transformation
     402 */
     403psPlaneTransform* psPlaneTransformCombine(
     404    psPlaneTransform *out,             ///< a transform to recycle, or NULL if one is to be created.
     405    const psPlaneTransform *trans1,    ///< first transform to combine
     406    const psPlaneTransform *trans2     ///< first transform to combine
     407);
     408
     409
     410/** takes two arrays containing matched coordinates and returns the
     411 *  best-fitting transformation.
     412 *
     413 *  psPlaneTransformFit takes two arrays containing matched coordinates (i.e.,
     414 *  coordinates in the source array correspond to the coordinates in the dest
     415 *  array) and returns the best-fitting transformation. The source and dest
     416 *  will contain psCoords. In the event that the number of coordinates in each
     417 *  is not identical, the function shall generate a warning, and extra
     418 *  coordinates in the longer of the two shall be ignored. The trans transform
     419 *  may not be NULL, since it specifies the desired order, polynomial type and
     420 *  any polynomial terms to mask. nRejIter rejection iterations shall be
     421 *  performed, wherein coordinates lying more than sigmaClip standard
     422 *  deviations from the fit shall be rejected.
     423 *
     424 *  @return bool        TRUE if successful, otherwise FALSE.
     425 */
     426bool psPlaneTransformFit(
     427    psPlaneTransform *trans,
     428    const psArray *source,
     429    const psArray *dest,
     430    int nRejIter,
     431    float sigmaClip
     432);
    350433
    351434/// @}

  • trunk/psLib/src/astronomy/psAstrometry.c

    r3559 r3598  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-03-30 02:21:14 $
     10 *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-03-31 23:01:46 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5151}
    5252
    53 /*****************************************************************************
    54 isProjectionLinear(): this is a private function which simply determines
    55 if the supplied psPlaneTransform transform is linear: if any of the
    56 cooefficients of order 2 are higher are non-zero, then it is not linear.
    57  *****************************************************************************/
    58 static psS32 isProjectionLinear(psPlaneTransform *transform)
    59 {
    60     PS_PTR_CHECK_NULL(transform, 0);
    61     PS_PTR_CHECK_NULL(transform->x, 0);
    62     PS_PTR_CHECK_NULL(transform->y, 0);
    63 
    64     for (psS32 i=0;i<(transform->x->nX);i++) {
    65         for (psS32 j=0;j<(transform->x->nY);j++) {
    66             if (transform->x->coeff[i][j] != 0.0) {
    67                 if (!(((i == 0) && (j == 0)) ||
    68                         ((i == 0) && (j == 1)) ||
    69                         ((i == 1) && (j == 0)))) {
    70                     return(0);
    71                 }
    72             }
    73         }
    74     }
    75 
    76     for (psS32 i=0;i<(transform->y->nX);i++) {
    77         for (psS32 j=0;j<(transform->y->nY);j++) {
    78             if (transform->y->coeff[i][j] != 0.0) {
    79                 if (!(((i == 0) && (j == 0)) ||
    80                         ((i == 0) && (j == 1)) ||
    81                         ((i == 1) && (j == 0)))) {
    82                     return(0);
    83                 }
    84             }
    85         }
    86     }
    87 
    88     return(1);
    89 }
    90 
    91 /*****************************************************************************
    92 invertPlaneTransform(transform): : this is a private function which
    93 simply inverts the supplied psPlaneTransform transform.  It assumes that
    94 "transform" is linear.
    95  
    96 This program assumes that the inverse of the following linear equations:
    97         X2 = A + (B * X1) + (C * Y1);
    98         Y2 = D + (E * X1) + (F * Y1);
    99 is
    100         Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B));
    101         X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
    102 or
    103  X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) +
    104       (X2 * -((F/C) / (E - ((F*B)/C)))) +
    105              (Y2 * (1.0 / (E - ((F*B)/C))));
    106         Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) +
    107              (X2 * -((E/B) / (F - ((C*E)/B)))) +
    108              (Y2 * (1.0 / (F - ((C*E)/B))));
    109  
    110 XXX: Since thre is now a general psPlaneTransformInvert() function, we
    111 should rename this.
    112  
    113  *****************************************************************************/
    114 static psPlaneTransform *invertPlaneTransform(psPlaneTransform *transform)
    115 {
    116     PS_PTR_CHECK_NULL(transform, 0);
    117     PS_PTR_CHECK_NULL(transform->x, 0);
    118     PS_PTR_CHECK_NULL(transform->y, 0);
    119 
    120     psF64 A = 0.0;
    121     psF64 B = 0.0;
    122     psF64 C = 0.0;
    123     psF64 D = 0.0;
    124     psF64 E = 0.0;
    125     psF64 F = 0.0;
    126 
    127     // XXX: Test this for correctness.
    128     A = transform->x->coeff[0][0];
    129     if (transform->x->nX >= 2) {
    130         B = transform->x->coeff[1][0];
    131     }
    132     if (transform->x->nY >= 2) {
    133         C = transform->x->coeff[0][1];
    134     }
    135     D = transform->y->coeff[0][0];
    136     if (transform->y->nX >= 2) {
    137         E = transform->y->coeff[1][0];
    138     }
    139     if (transform->y->nY >= 2) {
    140         F = transform->y->coeff[0][1];
    141     }
    142 
    143     // XXX: Use the constructor here.
    144     psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
    145 
    146     out->x->coeff[0][0] = -D + ((F*A)/C) / (E - ((F*B)/C));
    147     out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C));
    148     out->x->coeff[0][1] =  1.0 / (E - ((F*B)/C));
    149     out->y->coeff[0][0] = -D + ((E*A)/B) / (F - ((C*E)/B));
    150     out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B));
    151     out->y->coeff[0][1] =  1.0 / (F - ((C*E)/B));
    152 
    153     return(out);
    154 }
    15553
    15654static void FPAFree(psFPA* fpa)
    … …  
    914812
    915813    // generate an error if cell->toTP is not linear.
    916     if (0 == isProjectionLinear(cell->toTP)) {
     814    if (0 == p_psIsProjectionLinear(cell->toTP)) {
    917815        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    918816                PS_ERRORTEXT_psAstrometry_NONLINEAR_TRANSFORM,
    … …  
    920818    }
    921819
    922     TPtoCell = invertPlaneTransform(cell->toTP);
     820    TPtoCell = p_psPlaneTransformLinearInvert(cell->toTP);
    923821    cellCoord = psPlaneTransformApply(cellCoord, TPtoCell, tpCoord);
    924822
    … …  
    930828
    931829
    932 /*****************************************************************************
    933 multiplyCoeffs(trans1, trans2): Takes two 2-D polynomials as input and
    934 multiplies them.  Basically, for each non-zero coeff in the trans1 coeff[][]
    935 array, you must multiply by all non-zero coeffs in trans2.
    936  
    937 XXX: Inefficient in that the out polynomial is allocated every time.
    938  *****************************************************************************/
    939 psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1,
    940                                  psDPolynomial2D *trans2)
    941 {
    942     psS32 orderX = (trans1->nX + trans2->nX) - 1;
    943     psS32 orderY = (trans1->nX + trans2->nX) - 1;
    944 
    945     psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
    946     for (psS32 i = 0 ; i < out->nX; i++) {
    947         for (psS32 j = 0 ; j < out->nY; j++) {
    948             out->coeff[i][j] = 0.0;
    949             out->mask[i][j] = 0;
    950         }
    951     }
    952 
    953     for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) {
    954         for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) {
    955             if (0.0 != trans1->coeff[t1x][t1y]) {
    956                 for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
    957                     for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
    958                         out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
    959                     }
    960                 }
    961             }
    962         }
    963     }
    964     return(out);
    965 }
    966 
    967 
    968 
    969 
    970 /*****************************************************************************
    971 psPlaneTransformCombine(out, trans1, trans2)
    972  
    973 XXX: Much room for optimization.  Currently, we call the polyMultiply
    974 routine far too many times.
    975  *****************************************************************************/
    976 psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out,
    977         const psPlaneTransform *trans1,
    978         const psPlaneTransform *trans2)
    979 {
    980     PS_PTR_CHECK_NULL(trans1, NULL);
    981     PS_PTR_CHECK_NULL(trans2, NULL);
    982 
    983     //
    984     // Determine the size of the new psPlaneTransform.
    985     //
    986     // PS_MAX(  Number of x terms in T2->x * number of x terms in T1->x,
    987     //          Number of y terms in T2->x * number of x terms in T1->y,
    988     psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX),
    989                             (trans2->x->nY * trans1->y->nX));
    990     psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY),
    991                             (trans2->x->nY * trans1->y->nY));
    992 
    993     psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX),
    994                             (trans2->y->nY * trans1->y->nX));
    995     psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY),
    996                             (trans2->y->nY * trans1->y->nY));
    997     psS32 orderX = PS_MAX(orderXnX, orderYnX);
    998     psS32 orderY = PS_MAX(orderXnY, orderYnY);
    999 
    1000     //
    1001     // Allocate the new psPlaneTransform, if necessary.
    1002     //
    1003     psPlaneTransform *myPT = NULL;
    1004     if (out == NULL) {
    1005         myPT = psPlaneTransformAlloc(orderX, orderY);
    1006     } else {
    1007         if ((out->x->nX == orderX) && (out->x->nY == orderY) &&
    1008                 (out->y->nX == orderX) && (out->y->nY == orderY)) {
    1009             myPT = out;
    1010         } else {
    1011             psFree(out);
    1012             myPT = psPlaneTransformAlloc(orderX, orderY);
    1013         }
    1014     }
    1015 
    1016     //
    1017     // Initialize the new psPlaneTransform, if necessary.
    1018     //
    1019     for (psS32 i = 0 ; i < orderX ; i++) {
    1020         for (psS32 j = 0 ; j < orderY ; j++) {
    1021             myPT->x->coeff[i][j] = 0.0;
    1022             myPT->x->mask[i][j] = 0;
    1023             myPT->y->coeff[i][j] = 0.0;
    1024             myPT->y->mask[i][j] = 0;
    1025         }
    1026     }
    1027 
    1028     //
    1029     // For each term (a * x^i * y^j) in trans2, we substitute the appropriate
    1030     // equation from trans1, and raise it to the appropriate power.  This is
    1031     // done via the multiplyDPoly2D().  The result is a polynomial (currPoly)
    1032     // and its coefficients are added into the myPT coeff matrix.
    1033     //
    1034     // XXX: This is horribly inefficient in that the trans1 polys are repeatedly
    1035     // multiplied against themselves.  This can easily be improved.
    1036     //
    1037     for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) {
    1038         for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) {
    1039             psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
    1040             currPoly->coeff[0][0] = 1.0;
    1041             currPoly->mask[0][0] = 0;
    1042             psDPolynomial2D *newPoly = NULL;
    1043 
    1044             if (trans2->x->mask[t2x][t2y] == 0) {
    1045 
    1046                 // Must raise trans1->y to the t2y-power.
    1047                 for (psS32 c = 0 ; c < t2y; c++) {
    1048                     newPoly = multiplyDPoly2D(currPoly, trans1->y);
    1049                     psFree(currPoly);
    1050                     currPoly = newPoly;
    1051                 }
    1052 
    1053                 // Must raise trans1->x to the t2x-power.
    1054                 for (psS32 c = 0 ; c < t2x; c++) {
    1055                     newPoly = multiplyDPoly2D(currPoly, trans1->x);
    1056                     psFree(currPoly);
    1057                     currPoly = newPoly;
    1058                 }
    1059 
    1060                 // Set the appropriate coeffs in myPT->x
    1061                 for (psS32 i = 0 ; i < currPoly->nX ; i++) {
    1062                     for (psS32 j = 0 ; j < currPoly->nY ; j++) {
    1063                         myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
    1064                     }
    1065                 }
    1066             }
    1067             psFree(currPoly);
    1068         }
    1069     }
    1070 
    1071 
    1072 
    1073     for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) {
    1074         for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) {
    1075             psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
    1076             currPoly->coeff[0][0] = 1.0;
    1077             currPoly->mask[0][0] = 0;
    1078             psDPolynomial2D *newPoly = NULL;
    1079 
    1080             if (trans2->y->mask[t2x][t2y] == 0) {
    1081 
    1082                 // Must raise trans1->y to the t2y-power.
    1083                 for (psS32 c = 0 ; c < t2y; c++) {
    1084                     newPoly = multiplyDPoly2D(currPoly, trans1->y);
    1085                     psFree(currPoly);
    1086                     currPoly = newPoly;
    1087                 }
    1088 
    1089                 // Must raise trans1->x to the t2x-power.
    1090                 for (psS32 c = 0 ; c < t2x; c++) {
    1091                     newPoly = multiplyDPoly2D(currPoly, trans1->x);
    1092                     psFree(currPoly);
    1093                     currPoly = newPoly;
    1094                 }
    1095 
    1096                 // Set the appropriate coeffs in myPT->x
    1097                 for (psS32 i = 0 ; i < currPoly->nX ; i++) {
    1098                     for (psS32 j = 0 ; j < currPoly->nY ; j++) {
    1099                         myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
    1100                     }
    1101                 }
    1102             }
    1103             psFree(currPoly);
    1104         }
    1105     }
    1106 
    1107     return(myPT);
    1108 }
    1109 
    1110 /*****************************************************************************
    1111 psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip)
    1112  
    1113 XXX: What about nRejIter?  Iterations?
    1114 XXX: Use static vectors for internal data.
    1115  *****************************************************************************/
    1116 bool psPlaneTransformFit(psPlaneTransform *trans,
    1117                          const psArray *source,
    1118                          const psArray *dest,
    1119                          int nRejIter,
    1120                          float sigmaClip)
    1121 {
    1122     PS_PTR_CHECK_NULL(trans, NULL);
    1123     PS_PTR_CHECK_NULL(source, NULL);
    1124     PS_PTR_CHECK_NULL(dest, NULL);
    1125 
    1126     psS32 numCoords = PS_MIN(source->n, dest->n);
    1127     // This is not really necessary because of above conditionals.
    1128     psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
    1129 
    1130     //
    1131     // Create fake polynomial to use in evaluation
    1132     //
    1133     psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
    1134     for (int i = 0; i < order; i++) {
    1135         for (int j = 0; j < order; j++) {
    1136             fakePoly->coeff[i][j] = 1.0;
    1137             fakePoly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
    1138         }
    1139     }
    1140 
    1141     //
    1142     // Initialize the matrix and vectors
    1143     //
    1144     psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
    1145     psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
    1146     psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
    1147     psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
    1148     for (psS32 i = 0; i < nCoeff; i++) {
    1149         for (psS32 j = 0; j < nCoeff; j++) {
    1150             matrix->data.F64[i][j] = 0.0;
    1151         }
    1152         xVector->data.F64[i] = 0.0;
    1153         yVector->data.F64[i] = 0.0;
    1154     }
    1155 
    1156     //
    1157     // Iterate over the grid points
    1158     //
    1159     for (psS32 g = 0; g < numCoords; g++) {
    1160         // Iterate over the polynomial coefficients, accumulating the matrix and vectors
    1161 
    1162         for (psS32 i = 0, ijIndex = 0; i < order; i++) {
    1163             for (psS32 j = 0; j < order - i; j++, ijIndex++) {
    1164                 fakePoly->mask[i][j] = 0;
    1165                 psF64 xIn = ((psPlane *) source->data[g])->x;
    1166                 psF64 yIn = ((psPlane *) source->data[g])->y;
    1167                 psF64 xOut = ((psPlane *) dest->data[g])->x;
    1168                 psF64 yOut = ((psPlane *) dest->data[g])->y;
    1169                 psF64 ijPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
    1170                 fakePoly->mask[i][j] = 1;
    1171 
    1172                 for (psS32 m = 0, mnIndex = 0; m < order; m++) {
    1173                     for (psS32 n = 0; n < order - m; n++, mnIndex++) {
    1174                         fakePoly->mask[m][n] = 0;
    1175                         psF64 mnPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
    1176                         fakePoly->mask[m][n] = 1;
    1177 
    1178                         matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
    1179                     }
    1180                 }
    1181 
    1182                 xVector->data.F64[ijIndex] += ijPoly * xOut;
    1183                 yVector->data.F64[ijIndex] += ijPoly * yOut;
    1184             }
    1185         }
    1186     }
    1187 
    1188     //
    1189     // Solution via LU Decomposition
    1190     //
    1191     psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
    1192     psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
    1193     psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
    1194     psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
    1195 
    1196     //
    1197     // XXX: Should check the output of the matrix routines and return false if bad.
    1198     //
    1199 
    1200     //
    1201     // Stuff coefficients into transformation
    1202     //
    1203     for (psS32 i = 0, ijIndex = 0; i < order; i++) {
    1204         for (psS32 j = 0; j < order - i; j++, ijIndex++) {
    1205             trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
    1206             trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
    1207         }
    1208     }
    1209 
    1210     psFree(fakePoly);
    1211     psFree(permutation);
    1212     psFree(luMatrix);
    1213     psFree(xSolution);
    1214     psFree(ySolution);
    1215     psFree(matrix);
    1216     psFree(xVector);
    1217     psFree(yVector);
    1218 
    1219     return(true);
    1220 }
    1221 
    1222 
    1223 /*****************************************************************************
    1224 psPlaneTransformInvert(out, in, region, nSamples)
    1225  
    1226 // XXX: Use static data structures.
    1227  *****************************************************************************/
    1228 psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
    1229         const psPlaneTransform *in,
    1230         psRegion *region,
    1231         int nSamples)
    1232 {
    1233     PS_PTR_CHECK_NULL(in, NULL);
    1234     //
    1235     // If the transform is linear, then invert it exactly and return.
    1236     //
    1237     if (isProjectionLinear((psPlaneTransform *) in)) {
    1238         printf("COOL: is linear\n");
    1239         return(invertPlaneTransform((psPlaneTransform *) in));
    1240     }
    1241     PS_PTR_CHECK_NULL(region, NULL);
    1242     PS_INT_COMPARE(1, nSamples, NULL);
    1243 
    1244     // Ensure that the input transformation is symmetrical.
    1245     if ((in->x->nX != in->x->nY) ||
    1246             (in->y->nX != in->y->nY) ||
    1247             (in->x->nX != in->y->nX)) {
    1248         psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
    1249     }
    1250     psS32 order = PS_MAX(in->x->nX, in->x->nY);
    1251 
    1252     psPlaneTransform *myPT = NULL;
    1253     psPlane *inCoord = psPlaneAlloc();
    1254     psPlane *outCoord = psPlaneAlloc();
    1255 
    1256     //
    1257     // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
    1258     //
    1259     if (out == NULL) {
    1260         myPT = psPlaneTransformAlloc(order, order);
    1261     } else {
    1262         if ((out->x->nX == order) && (out->x->nY == order) &&
    1263                 (out->y->nX == order) && (out->y->nY == order)) {
    1264             myPT = out;
    1265         } else {
    1266             psFree(out);
    1267             myPT = psPlaneTransformAlloc(order, order);
    1268         }
    1269     }
    1270 
    1271     //
    1272     // Copy the input transform to myPT.
    1273     //
    1274     for (psS32 i = 0 ; i < in->x->nX ; i++) {
    1275         for (psS32 j = 0 ; j < in->x->nY ; j++) {
    1276             myPT->x->coeff[i][j] = in->x->coeff[i][j];
    1277         }
    1278     }
    1279     for (psS32 i = 0 ; i < in->y->nX ; i++) {
    1280         for (psS32 j = 0 ; j < in->y->nY ; j++) {
    1281             myPT->y->coeff[i][j] = in->y->coeff[i][j];
    1282         }
    1283     }
    1284 
    1285     //
    1286     // Create a grid of xin,yin --> xout,yout
    1287     //
    1288     psArray *inData = psArrayAlloc(nSamples * nSamples);
    1289     psArray *outData = psArrayAlloc(nSamples * nSamples);
    1290     for (psS32 i = 0 ; i < inData->n; i++) {
    1291         inData->data[i] = (psPtr *) psPlaneAlloc();
    1292         outData->data[i] = (psPtr *) psPlaneAlloc();
    1293     }
    1294 
    1295     //
    1296     // Initialize the grid.
    1297     //
    1298     psS32 cnt = 0;
    1299     for (int yint = 0; yint < nSamples; yint++) {
    1300         inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
    1301         for (int xint = 0; xint < nSamples; xint++) {
    1302             inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
    1303             (void)psPlaneTransformApply(outCoord, in, inCoord);
    1304 
    1305             ((psPlane *) outData->data[cnt])->x = inCoord->x;
    1306             ((psPlane *) outData->data[cnt])->y = inCoord->y;
    1307             ((psPlane *) inData->data[cnt])->x = outCoord->x;
    1308             ((psPlane *) inData->data[cnt])->y = outCoord->y;
    1309 
    1310             cnt++;
    1311         }
    1312     }
    1313     bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0);
    1314 
    1315     psFree(inCoord);
    1316     psFree(outCoord);
    1317     psFree(inData);
    1318     psFree(outData);
    1319 
    1320     if (rc == true) {
    1321         return(myPT);
    1322     }
    1323 
    1324     // XXX: Generate an error message, or warning message.
    1325     return(NULL);
    1326 }
     830

  • trunk/psLib/src/astronomy/psAstrometry.h

    r3587 r3598  
    88*  @author GLG, MHPCC
    99*
    10 *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2005-03-31 02:33:28 $
     10*  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2005-03-31 23:01:46 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    535535
    536536
    537 // XXX: These functions don't belong here.  Will migrate to psCoords.c later.
    538 // XXX: Doxygenate.
    539 psPlaneTransform *psPlaneTransformInvert(
    540     psPlaneTransform *out,
    541     const psPlaneTransform *in,
    542     psRegion *region,
    543     int nSamples
    544 );
    545 
    546 // XXX: Doxygenate.
    547 psPlaneTransform *psPlaneTransformCombine(
    548     psPlaneTransform *out,
    549     const psPlaneTransform *trans1,
    550     const psPlaneTransform *trans2
    551 );
    552 
    553 // XXX: Doxygenate.
    554 bool psPlaneTransformFit(
    555     psPlaneTransform *trans,
    556     const psArray *source,
    557     const psArray *dest,
    558     int nRejIter,
    559     float sigmaClip
    560 );
    561 
    562537#endif

  • trunk/psLib/src/astronomy/psCoord.c

    r3540 r3598  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-29 19:41:56 $
     12*  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-31 23:01:46 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2626#include "psLogMsg.h"
    2727#include "psAstronomyErrors.h"
     28#include "psAstrometry.h"
     29#include "psMatrix.h"
    2830#include <math.h>
    2931#include <float.h>
    … …  
    5961/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
    6062/*****************************************************************************/
     63/*****************************************************************************
     64multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and
     65multiplies them.  Basically, for each non-zero coeff in the trans1 coeff[][]
     66array, you must multiply by all non-zero coeffs in trans2.
     67 
     68XXX: Inefficient in that the out polynomial is allocated every time.
     69 *****************************************************************************/
     70
     71psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1,
     72                                 psDPolynomial2D *trans2)
     73{
     74    psS32 orderX = (trans1->nX + trans2->nX) - 1;
     75    psS32 orderY = (trans1->nX + trans2->nX) - 1;
     76
     77    psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
     78    for (psS32 i = 0 ; i < out->nX; i++) {
     79        for (psS32 j = 0 ; j < out->nY; j++) {
     80            out->coeff[i][j] = 0.0;
     81            out->mask[i][j] = 0;
     82        }
     83    }
     84
     85    for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) {
     86        for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) {
     87            if (0.0 != trans1->coeff[t1x][t1y]) {
     88                for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
     89                    for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
     90                        out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
     91                    }
     92                }
     93            }
     94        }
     95    }
     96    return(out);
     97}
    6198
    6299/*****************************************************************************/
    … …  
    66103{
    67104    // There are non dynamic allocated items
     105}
     106
     107/*****************************************************************************
     108p_psPlaneTransformLinearInvert(transform): : this is a private function which
     109simply inverts the supplied psPlaneTransform transform.  It assumes that
     110"transform" is linear.
     111 
     112This program assumes that the inverse of the following linear equations:
     113        X2 = A + (B * X1) + (C * Y1);
     114        Y2 = D + (E * X1) + (F * Y1);
     115is
     116        Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B));
     117        X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
     118or
     119 X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) +
     120      (X2 * -((F/C) / (E - ((F*B)/C)))) +
     121             (Y2 * (1.0 / (E - ((F*B)/C))));
     122        Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) +
     123             (X2 * -((E/B) / (F - ((C*E)/B)))) +
     124             (Y2 * (1.0 / (F - ((C*E)/B))));
     125 
     126XXX: Since thre is now a general psPlaneTransformInvert() function, we
     127should rename this.
     128 
     129 *****************************************************************************/
     130psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
     131{
     132    PS_PTR_CHECK_NULL(transform, 0);
     133    PS_PTR_CHECK_NULL(transform->x, 0);
     134    PS_PTR_CHECK_NULL(transform->y, 0);
     135
     136    psF64 A = 0.0;
     137    psF64 B = 0.0;
     138    psF64 C = 0.0;
     139    psF64 D = 0.0;
     140    psF64 E = 0.0;
     141    psF64 F = 0.0;
     142
     143    // XXX: Test this for correctness.
     144    A = transform->x->coeff[0][0];
     145    if (transform->x->nX >= 2) {
     146        B = transform->x->coeff[1][0];
     147    }
     148    if (transform->x->nY >= 2) {
     149        C = transform->x->coeff[0][1];
     150    }
     151    D = transform->y->coeff[0][0];
     152    if (transform->y->nX >= 2) {
     153        E = transform->y->coeff[1][0];
     154    }
     155    if (transform->y->nY >= 2) {
     156        F = transform->y->coeff[0][1];
     157    }
     158
     159    // XXX: Use the constructor here.
     160    psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
     161
     162    out->x->coeff[0][0] = -D + ((F*A)/C) / (E - ((F*B)/C));
     163    out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C));
     164    out->x->coeff[0][1] =  1.0 / (E - ((F*B)/C));
     165    out->y->coeff[0][0] = -D + ((E*A)/B) / (F - ((C*E)/B));
     166    out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B));
     167    out->y->coeff[0][1] =  1.0 / (F - ((C*E)/B));
     168
     169    return(out);
     170}
     171
     172/*****************************************************************************
     173p_psIsProjectionLinear(): this is a private function which simply determines
     174if the supplied psPlaneTransform transform is linear: if any of the
     175cooefficients of order 2 are higher are non-zero, then it is not linear.
     176 *****************************************************************************/
     177psS32 p_psIsProjectionLinear(psPlaneTransform *transform)
     178{
     179    PS_PTR_CHECK_NULL(transform, 0);
     180    PS_PTR_CHECK_NULL(transform->x, 0);
     181    PS_PTR_CHECK_NULL(transform->y, 0);
     182
     183    for (psS32 i=0;i<(transform->x->nX);i++) {
     184        for (psS32 j=0;j<(transform->x->nY);j++) {
     185            if (transform->x->coeff[i][j] != 0.0) {
     186                if (!(((i == 0) && (j == 0)) ||
     187                        ((i == 0) && (j == 1)) ||
     188                        ((i == 1) && (j == 0)))) {
     189                    return(0);
     190                }
     191            }
     192        }
     193    }
     194
     195    for (psS32 i=0;i<(transform->y->nX);i++) {
     196        for (psS32 j=0;j<(transform->y->nY);j++) {
     197            if (transform->y->coeff[i][j] != 0.0) {
     198                if (!(((i == 0) && (j == 0)) ||
     199                        ((i == 0) && (j == 1)) ||
     200                        ((i == 1) && (j == 0)))) {
     201                    return(0);
     202                }
     203            }
     204        }
     205    }
     206
     207    return(1);
    68208}
    69209
    … …  
    721861}
    722862
     863
     864
     865/*****************************************************************************
     866psPlaneTransformCombine(out, trans1, trans2)
     867 
     868XXX: Much room for optimization.  Currently, we call the polyMultiply
     869routine far too many times.
     870 *****************************************************************************/
     871psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out,
     872        const psPlaneTransform *trans1,
     873        const psPlaneTransform *trans2)
     874{
     875    PS_PTR_CHECK_NULL(trans1, NULL);
     876    PS_PTR_CHECK_NULL(trans2, NULL);
     877
     878    //
     879    // Determine the size of the new psPlaneTransform.
     880    //
     881    // PS_MAX(  Number of x terms in T2->x * number of x terms in T1->x,
     882    //          Number of y terms in T2->x * number of x terms in T1->y,
     883    psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX),
     884                            (trans2->x->nY * trans1->y->nX));
     885    psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY),
     886                            (trans2->x->nY * trans1->y->nY));
     887
     888    psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX),
     889                            (trans2->y->nY * trans1->y->nX));
     890    psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY),
     891                            (trans2->y->nY * trans1->y->nY));
     892    psS32 orderX = PS_MAX(orderXnX, orderYnX);
     893    psS32 orderY = PS_MAX(orderXnY, orderYnY);
     894
     895    //
     896    // Allocate the new psPlaneTransform, if necessary.
     897    //
     898    psPlaneTransform *myPT = NULL;
     899    if (out == NULL) {
     900        myPT = psPlaneTransformAlloc(orderX, orderY);
     901    } else {
     902        if ((out->x->nX == orderX) && (out->x->nY == orderY) &&
     903                (out->y->nX == orderX) && (out->y->nY == orderY)) {
     904            myPT = out;
     905        } else {
     906            psFree(out);
     907            myPT = psPlaneTransformAlloc(orderX, orderY);
     908        }
     909    }
     910
     911    //
     912    // Initialize the new psPlaneTransform, if necessary.
     913    //
     914    for (psS32 i = 0 ; i < orderX ; i++) {
     915        for (psS32 j = 0 ; j < orderY ; j++) {
     916            myPT->x->coeff[i][j] = 0.0;
     917            myPT->x->mask[i][j] = 0;
     918            myPT->y->coeff[i][j] = 0.0;
     919            myPT->y->mask[i][j] = 0;
     920        }
     921    }
     922
     923    //
     924    // For each term (a * x^i * y^j) in trans2, we substitute the appropriate
     925    // equation from trans1, and raise it to the appropriate power.  This is
     926    // done via the multiplyDPoly2D().  The result is a polynomial (currPoly)
     927    // and its coefficients are added into the myPT coeff matrix.
     928    //
     929    // XXX: This is horribly inefficient in that the trans1 polys are repeatedly
     930    // multiplied against themselves.  This can easily be improved.
     931    //
     932    for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) {
     933        for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) {
     934            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
     935            currPoly->coeff[0][0] = 1.0;
     936            currPoly->mask[0][0] = 0;
     937            psDPolynomial2D *newPoly = NULL;
     938
     939            if (trans2->x->mask[t2x][t2y] == 0) {
     940
     941                // Must raise trans1->y to the t2y-power.
     942                for (psS32 c = 0 ; c < t2y; c++) {
     943                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
     944                    psFree(currPoly);
     945                    currPoly = newPoly;
     946                }
     947
     948                // Must raise trans1->x to the t2x-power.
     949                for (psS32 c = 0 ; c < t2x; c++) {
     950                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
     951                    psFree(currPoly);
     952                    currPoly = newPoly;
     953                }
     954
     955                // Set the appropriate coeffs in myPT->x
     956                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
     957                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
     958                        myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
     959                    }
     960                }
     961            }
     962            psFree(currPoly);
     963        }
     964    }
     965
     966
     967
     968    for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) {
     969        for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) {
     970            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
     971            currPoly->coeff[0][0] = 1.0;
     972            currPoly->mask[0][0] = 0;
     973            psDPolynomial2D *newPoly = NULL;
     974
     975            if (trans2->y->mask[t2x][t2y] == 0) {
     976
     977                // Must raise trans1->y to the t2y-power.
     978                for (psS32 c = 0 ; c < t2y; c++) {
     979                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
     980                    psFree(currPoly);
     981                    currPoly = newPoly;
     982                }
     983
     984                // Must raise trans1->x to the t2x-power.
     985                for (psS32 c = 0 ; c < t2x; c++) {
     986                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
     987                    psFree(currPoly);
     988                    currPoly = newPoly;
     989                }
     990
     991                // Set the appropriate coeffs in myPT->x
     992                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
     993                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
     994                        myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
     995                    }
     996                }
     997            }
     998            psFree(currPoly);
     999        }
     1000    }
     1001
     1002    return(myPT);
     1003}
     1004
     1005/*****************************************************************************
     1006psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip)
     1007 
     1008XXX: What about nRejIter?  Iterations?
     1009XXX: Use static vectors for internal data.
     1010 *****************************************************************************/
     1011bool psPlaneTransformFit(psPlaneTransform *trans,
     1012                         const psArray *source,
     1013                         const psArray *dest,
     1014                         int nRejIter,
     1015                         float sigmaClip)
     1016{
     1017    PS_PTR_CHECK_NULL(trans, NULL);
     1018    PS_PTR_CHECK_NULL(source, NULL);
     1019    PS_PTR_CHECK_NULL(dest, NULL);
     1020
     1021    psS32 numCoords = PS_MIN(source->n, dest->n);
     1022    // This is not really necessary because of above conditionals.
     1023    psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
     1024
     1025    //
     1026    // Create fake polynomial to use in evaluation
     1027    //
     1028    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
     1029    for (int i = 0; i < order; i++) {
     1030        for (int j = 0; j < order; j++) {
     1031            fakePoly->coeff[i][j] = 1.0;
     1032            fakePoly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
     1033        }
     1034    }
     1035
     1036    //
     1037    // Initialize the matrix and vectors
     1038    //
     1039    psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
     1040    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
     1041    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
     1042    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
     1043    for (psS32 i = 0; i < nCoeff; i++) {
     1044        for (psS32 j = 0; j < nCoeff; j++) {
     1045            matrix->data.F64[i][j] = 0.0;
     1046        }
     1047        xVector->data.F64[i] = 0.0;
     1048        yVector->data.F64[i] = 0.0;
     1049    }
     1050
     1051    //
     1052    // Iterate over the grid points
     1053    //
     1054    for (psS32 g = 0; g < numCoords; g++) {
     1055        // Iterate over the polynomial coefficients, accumulating the matrix and vectors
     1056
     1057        for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1058            for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1059                fakePoly->mask[i][j] = 0;
     1060                psF64 xIn = ((psPlane *) source->data[g])->x;
     1061                psF64 yIn = ((psPlane *) source->data[g])->y;
     1062                psF64 xOut = ((psPlane *) dest->data[g])->x;
     1063                psF64 yOut = ((psPlane *) dest->data[g])->y;
     1064                psF64 ijPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
     1065                fakePoly->mask[i][j] = 1;
     1066
     1067                for (psS32 m = 0, mnIndex = 0; m < order; m++) {
     1068                    for (psS32 n = 0; n < order - m; n++, mnIndex++) {
     1069                        fakePoly->mask[m][n] = 0;
     1070                        psF64 mnPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
     1071                        fakePoly->mask[m][n] = 1;
     1072
     1073                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
     1074                    }
     1075                }
     1076
     1077                xVector->data.F64[ijIndex] += ijPoly * xOut;
     1078                yVector->data.F64[ijIndex] += ijPoly * yOut;
     1079            }
     1080        }
     1081    }
     1082
     1083    //
     1084    // Solution via LU Decomposition
     1085    //
     1086    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     1087    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
     1088    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
     1089    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
     1090
     1091    //
     1092    // XXX: Should check the output of the matrix routines and return false if bad.
     1093    //
     1094
     1095    //
     1096    // Stuff coefficients into transformation
     1097    //
     1098    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1099        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1100            trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
     1101            trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
     1102        }
     1103    }
     1104
     1105    psFree(fakePoly);
     1106    psFree(permutation);
     1107    psFree(luMatrix);
     1108    psFree(xSolution);
     1109    psFree(ySolution);
     1110    psFree(matrix);
     1111    psFree(xVector);
     1112    psFree(yVector);
     1113
     1114    return(true);
     1115}
     1116
     1117
     1118/*****************************************************************************
     1119psPlaneTransformInvert(out, in, region, nSamples)
     1120 
     1121// XXX: Use static data structures.
     1122 *****************************************************************************/
     1123psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
     1124        const psPlaneTransform *in,
     1125        psRegion *region,
     1126        int nSamples)
     1127{
     1128    PS_PTR_CHECK_NULL(in, NULL);
     1129    //
     1130    // If the transform is linear, then invert it exactly and return.
     1131    //
     1132    if (p_psIsProjectionLinear((psPlaneTransform *) in)) {
     1133        printf("COOL: is linear\n");
     1134        return(p_psPlaneTransformLinearInvert((psPlaneTransform *) in));
     1135    }
     1136    PS_PTR_CHECK_NULL(region, NULL);
     1137    PS_INT_COMPARE(1, nSamples, NULL);
     1138
     1139    // Ensure that the input transformation is symmetrical.
     1140    if ((in->x->nX != in->x->nY) ||
     1141            (in->y->nX != in->y->nY) ||
     1142            (in->x->nX != in->y->nX)) {
     1143        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
     1144    }
     1145    psS32 order = PS_MAX(in->x->nX, in->x->nY);
     1146
     1147    psPlaneTransform *myPT = NULL;
     1148    psPlane *inCoord = psPlaneAlloc();
     1149    psPlane *outCoord = psPlaneAlloc();
     1150
     1151    //
     1152    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
     1153    //
     1154    if (out == NULL) {
     1155        myPT = psPlaneTransformAlloc(order, order);
     1156    } else {
     1157        if ((out->x->nX == order) && (out->x->nY == order) &&
     1158                (out->y->nX == order) && (out->y->nY == order)) {
     1159            myPT = out;
     1160        } else {
     1161            psFree(out);
     1162            myPT = psPlaneTransformAlloc(order, order);
     1163        }
     1164    }
     1165
     1166    //
     1167    // Copy the input transform to myPT.
     1168    //
     1169    for (psS32 i = 0 ; i < in->x->nX ; i++) {
     1170        for (psS32 j = 0 ; j < in->x->nY ; j++) {
     1171            myPT->x->coeff[i][j] = in->x->coeff[i][j];
     1172        }
     1173    }
     1174    for (psS32 i = 0 ; i < in->y->nX ; i++) {
     1175        for (psS32 j = 0 ; j < in->y->nY ; j++) {
     1176            myPT->y->coeff[i][j] = in->y->coeff[i][j];
     1177        }
     1178    }
     1179
     1180    //
     1181    // Create a grid of xin,yin --> xout,yout
     1182    //
     1183    psArray *inData = psArrayAlloc(nSamples * nSamples);
     1184    psArray *outData = psArrayAlloc(nSamples * nSamples);
     1185    for (psS32 i = 0 ; i < inData->n; i++) {
     1186        inData->data[i] = (psPtr *) psPlaneAlloc();
     1187        outData->data[i] = (psPtr *) psPlaneAlloc();
     1188    }
     1189
     1190    //
     1191    // Initialize the grid.
     1192    //
     1193    psS32 cnt = 0;
     1194    for (int yint = 0; yint < nSamples; yint++) {
     1195        inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
     1196        for (int xint = 0; xint < nSamples; xint++) {
     1197            inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
     1198            (void)psPlaneTransformApply(outCoord, in, inCoord);
     1199
     1200            ((psPlane *) outData->data[cnt])->x = inCoord->x;
     1201            ((psPlane *) outData->data[cnt])->y = inCoord->y;
     1202            ((psPlane *) inData->data[cnt])->x = outCoord->x;
     1203            ((psPlane *) inData->data[cnt])->y = outCoord->y;
     1204
     1205            cnt++;
     1206        }
     1207    }
     1208    bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0);
     1209
     1210    psFree(inCoord);
     1211    psFree(outCoord);
     1212    psFree(inData);
     1213    psFree(outData);
     1214
     1215    if (rc == true) {
     1216        return(myPT);
     1217    }
     1218
     1219    // XXX: Generate an error message, or warning message.
     1220    return(NULL);
     1221}

  • trunk/psLib/src/astronomy/psCoord.h

    r3540 r3598  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-29 19:41:56 $
     12*  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-31 23:01:46 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    344344);
    345345
    346 psSphere *psSpherePrecess(psSphere *coords,
    347                           const psTime *fromTime,
    348                           const psTime *toTime);
    349 
     346/** Generates the complete spherical rotation to account for precession
     347 *  between two times.  The equinoxes shall be Julian equinoxes.
     348 *
     349 *  @return psSphere* the resulting spherical rotation
     350 */
     351psSphere* psSpherePrecess(
     352    psSphere *coords,                  ///< coordinates (modified in-place)
     353    const psTime *fromTime,            ///< equinox of coords input
     354    const psTime *toTime               ///< equinox of coords output
     355);
     356
     357// XXX: Doxygenate.
     358psPlaneTransform *p_psPlaneTransformLinearInvert(
     359    psPlaneTransform *transform
     360);
     361
     362// XXX: Doxygenate
     363psS32 p_psIsProjectionLinear(
     364    psPlaneTransform *transform
     365);
     366
     367/** inverts a given transformation.
     368 *
     369 *  It may assume that the input transformation is one-to-one, and that the
     370 *  inverse transformation may be specified through using polynomials of the
     371 *  same type and order as the forward transformation. In the event that the
     372 *  input transformation is linear, an exact solution may be calculated;
     373 *  otherwise nSamples samples in each axis, covering the region specified by
     374 *  region shall be used as a grid to fit the best inverse transformation. The
     375 *  function shall return NULL if it was unable to generate the inverse
     376 *  transformation; otherwise it shall return the inverse transformation. In
     377 *  the event that out is NULL, a new psPlaneTransform shall be allocated and
     378 *  returned.
     379 *
     380 *  @return psPlaneTransform*  the resulting inverted transform
     381 */
     382psPlaneTransform* psPlaneTransformInvert(
     383    psPlaneTransform *out,             ///< a transform to recycle, or NULL if one is to be created.
     384    const psPlaneTransform *in,        ///< transform to invert
     385    psRegion *region,                  ///< region to fit for non-linear transform inversion
     386    int nSamples                       ///< number of samples in each axis for fit
     387);
     388
     389/** Creates a single transformation that has the effect of performing trans1
     390 *  followed by trans2.
     391 *
     392 *  psPlaneTransformCombine takes two transformations (trans1 and trans2) and
     393 *  returns a single transformation that has the effect of performing trans1
     394 *  followed by trans2. In the event that the input transformation is linear,
     395 *  an exact solution may be calculated; otherwise nSamples samples in each
     396 *  axis, covering the region specified by region shall be used as a grid to
     397 *  fit the best inverse transformation. The function shall return NULL if it
     398 *  was unable to generate the transformation; otherwise it shall return the
     399 *  transformation.
     400 *
     401 *  @return psPlaneTransform*    resulting transformation
     402 */
     403psPlaneTransform* psPlaneTransformCombine(
     404    psPlaneTransform *out,             ///< a transform to recycle, or NULL if one is to be created.
     405    const psPlaneTransform *trans1,    ///< first transform to combine
     406    const psPlaneTransform *trans2     ///< first transform to combine
     407);
     408
     409
     410/** takes two arrays containing matched coordinates and returns the
     411 *  best-fitting transformation.
     412 *
     413 *  psPlaneTransformFit takes two arrays containing matched coordinates (i.e.,
     414 *  coordinates in the source array correspond to the coordinates in the dest
     415 *  array) and returns the best-fitting transformation. The source and dest
     416 *  will contain psCoords. In the event that the number of coordinates in each
     417 *  is not identical, the function shall generate a warning, and extra
     418 *  coordinates in the longer of the two shall be ignored. The trans transform
     419 *  may not be NULL, since it specifies the desired order, polynomial type and
     420 *  any polynomial terms to mask. nRejIter rejection iterations shall be
     421 *  performed, wherein coordinates lying more than sigmaClip standard
     422 *  deviations from the fit shall be rejected.
     423 *
     424 *  @return bool        TRUE if successful, otherwise FALSE.
     425 */
     426bool psPlaneTransformFit(
     427    psPlaneTransform *trans,
     428    const psArray *source,
     429    const psArray *dest,
     430    int nRejIter,
     431    float sigmaClip
     432);
    350433
    351434/// @}

  • trunk/psLib/src/astronomy/psDB.h

    r3591 r3598  
    1010 *  @author Joshua Hoblitt
    1111 *
    12  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-03-31 02:54:28 $
     12 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-03-31 23:01:46 $
    1414 *
    1515 *  Copyright 2005 Joshua Hoblitt, University of Hawaii
    … …  
    2121#ifdef BUILD_PSDB
    2222
    23 #include <mysql.h>
    24 
    2523#include "psType.h"
    2624#include "psMetadata.h"
    … …  
    3634typedef struct
    3735{
    38     MYSQL *mysql;   ///< MySQL database handle
     36    void* mysql;   ///< MySQL database handle
    3937}
    4038psDB;

  • trunk/psLib/src/dataIO/psDB.h

    r3591 r3598  
    1010 *  @author Joshua Hoblitt
    1111 *
    12  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-03-31 02:54:28 $
     12 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-03-31 23:01:46 $
    1414 *
    1515 *  Copyright 2005 Joshua Hoblitt, University of Hawaii
    … …  
    2121#ifdef BUILD_PSDB
    2222
    23 #include <mysql.h>
    24 
    2523#include "psType.h"
    2624#include "psMetadata.h"
    … …  
    3634typedef struct
    3735{
    38     MYSQL *mysql;   ///< MySQL database handle
     36    void* mysql;   ///< MySQL database handle
    3937}
    4038psDB;

  • trunk/psLib/src/dataManip/psConstants.h

    r3575 r3598  
    1 /** @file  psComments.h
     1/** @file  psConstants.h
    22 *
    33 *  This file will hold definitions of various constants as well as common
    … …  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-03-30 23:44:05 $
     8 *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-03-31 23:01:46 $
    1010 *
    1111 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1313 *  XXX: Add parenthesis around all arguments so that these macros can be
    1414 *       called with complex expressions.
    15  *
    16  *
    1715 *
    1816 */

  • trunk/psLib/src/dataManip/psFunctions.h

    r3264 r3598  
    11/** @file psFunctions.h
    2 *  \brief Standard Mathematical Functions.
    3 *  \ingroup Stats
    4 *
    5 *  This file will hold the prototypes for procedures which allocate, free,
    6 *  and evaluate various polynomials.  Those polynomial structures are also
    7 *  defined here.
    8 *
    9 *  @ingroup Stats
    10 *
    11 *  @author Someone at IfA
    12 *  @author GLG, MHPCC
    13 *
    14 *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2005-02-17 19:26:23 $
    16 *
    17 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    18 */
     2 *  @brief Standard Mathematical Functions.
     3 *  @ingroup Stats
     4 *
     5 *  This file will hold the prototypes for procedures which allocate, free,
     6 *  and evaluate various polynomials.  Those polynomial structures are also
     7 *  defined here.
     8 *
     9 *  @ingroup Stats
     10 *
     11 *  @author Someone at IfA
     12 *  @author GLG, MHPCC
     13 *
     14 *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2005-03-31 23:01:46 $
     16 *
     17 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     18 */
    1919
    2020#if !defined(PS_FUNCTIONS_H)

  • trunk/psLib/src/db/psDB.h

    r3591 r3598  
    1010 *  @author Joshua Hoblitt
    1111 *
    12  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-03-31 02:54:28 $
     12 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-03-31 23:01:46 $
    1414 *
    1515 *  Copyright 2005 Joshua Hoblitt, University of Hawaii
    … …  
    2121#ifdef BUILD_PSDB
    2222
    23 #include <mysql.h>
    24 
    2523#include "psType.h"
    2624#include "psMetadata.h"
    … …  
    3634typedef struct
    3735{
    38     MYSQL *mysql;   ///< MySQL database handle
     36    void* mysql;   ///< MySQL database handle
    3937}
    4038psDB;

  • trunk/psLib/src/fileUtils/psDB.h

    r3591 r3598  
    1010 *  @author Joshua Hoblitt
    1111 *
    12  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-03-31 02:54:28 $
     12 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-03-31 23:01:46 $
    1414 *
    1515 *  Copyright 2005 Joshua Hoblitt, University of Hawaii
    … …  
    2121#ifdef BUILD_PSDB
    2222
    23 #include <mysql.h>
    24 
    2523#include "psType.h"
    2624#include "psMetadata.h"
    … …  
    3634typedef struct
    3735{
    38     MYSQL *mysql;   ///< MySQL database handle
     36    void* mysql;   ///< MySQL database handle
    3937}
    4038psDB;

  • trunk/psLib/src/mainpage.dox

    r3584 r3598  
     1/** @file mainpage.dox
     2 *  @brief Main page for the Doxygen documentation
     3 *
     4 *  @author Robert Daniel DeSonia, MHPCC
     5 *
     6 *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2005-03-31 23:01:46 $
     8 *
     9 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     10 */
     11
    112/** @mainpage psLib Image Processing Library
    2 
    313
    414@section intro Introduction

  • trunk/psLib/src/math/psConstants.h

    r3575 r3598  
    1 /** @file  psComments.h
     1/** @file  psConstants.h
    22 *
    33 *  This file will hold definitions of various constants as well as common
    … …  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-03-30 23:44:05 $
     8 *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-03-31 23:01:46 $
    1010 *
    1111 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1313 *  XXX: Add parenthesis around all arguments so that these macros can be
    1414 *       called with complex expressions.
    15  *
    16  *
    1715 *
    1816 */

  • trunk/psLib/src/math/psPolynomial.h

    r3264 r3598  
    11/** @file psFunctions.h
    2 *  \brief Standard Mathematical Functions.
    3 *  \ingroup Stats
    4 *
    5 *  This file will hold the prototypes for procedures which allocate, free,
    6 *  and evaluate various polynomials.  Those polynomial structures are also
    7 *  defined here.
    8 *
    9 *  @ingroup Stats
    10 *
    11 *  @author Someone at IfA
    12 *  @author GLG, MHPCC
    13 *
    14 *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2005-02-17 19:26:23 $
    16 *
    17 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    18 */
     2 *  @brief Standard Mathematical Functions.
     3 *  @ingroup Stats
     4 *
     5 *  This file will hold the prototypes for procedures which allocate, free,
     6 *  and evaluate various polynomials.  Those polynomial structures are also
     7 *  defined here.
     8 *
     9 *  @ingroup Stats
     10 *
     11 *  @author Someone at IfA
     12 *  @author GLG, MHPCC
     13 *
     14 *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2005-03-31 23:01:46 $
     16 *
     17 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     18 */
    1919
    2020#if !defined(PS_FUNCTIONS_H)

  • trunk/psLib/src/math/psSpline.h

    r3264 r3598  
    11/** @file psFunctions.h
    2 *  \brief Standard Mathematical Functions.
    3 *  \ingroup Stats
    4 *
    5 *  This file will hold the prototypes for procedures which allocate, free,
    6 *  and evaluate various polynomials.  Those polynomial structures are also
    7 *  defined here.
    8 *
    9 *  @ingroup Stats
    10 *
    11 *  @author Someone at IfA
    12 *  @author GLG, MHPCC
    13 *
    14 *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2005-02-17 19:26:23 $
    16 *
    17 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    18 */
     2 *  @brief Standard Mathematical Functions.
     3 *  @ingroup Stats
     4 *
     5 *  This file will hold the prototypes for procedures which allocate, free,
     6 *  and evaluate various polynomials.  Those polynomial structures are also
     7 *  defined here.
     8 *
     9 *  @ingroup Stats
     10 *
     11 *  @author Someone at IfA
     12 *  @author GLG, MHPCC
     13 *
     14 *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2005-03-31 23:01:46 $
     16 *
     17 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     18 */
    1919
    2020#if !defined(PS_FUNCTIONS_H)

  • trunk/psLib/src/psTest.h

    r3115 r3598  
     1/** @file  psTest.h
     2 *
     3 *  Testing infrastructure functions.
     4 *
     5 *  @author Robert DeSonia, MHPCC
     6 *
     7 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-03-31 23:01:46 $
     9 *
     10 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     11 *
     12 */
     13
    114
    215#ifndef PSTEST_H

  • trunk/psLib/src/pslib.h

    r3592 r3598  
    99*  @author Eric Van Alst, MHPCC
    1010*
    11 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-03-31 03:02:15 $
     11*  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2005-03-31 23:01:46 $
    1313*
    1414*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2323
    2424/******************************************************************************/
    25 #include "psDB.h"
    2625
    2726// // System Utilities
    … …  
    169168#include "psAstrometry.h"
    170169#include "psConstants.h"
     170#include "psDB.h"
    171171
    172172/// @}

  • trunk/psLib/swig/Makefile.am

    r3168 r3598  
    2323
    2424psLibModule/Makefile.PL: psLibModule
    25         cp $(srcdir)/Makefile.PL psLibModule/Makefile.PL
     25        cp -f $(srcdir)/Makefile.PL psLibModule/Makefile.PL
    2626
    2727psLibModule/setup.txt:psLibModule Makefile
    … …  
    5050
    5151uninstall-local: psLibModule/Makefile
    52         cd PsLibModule && make clean
     52        cd psLibModule && make clean

  • trunk/psLib/test/astronomy/tst_psDB.c

    r3569 r3598  
    99 *  @author Aaron Culliney, MHPCC
    1010 *
    11  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-30 23:22:39 $
     11 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-31 23:01:46 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 */
    17 
    18 #define PS_ALLOW_MALLOC // for MySQL internals ...
    1917
    2018#include <unistd.h>
    … …  
    2624#include <string.h>
    2725
    28 #include "psDB.h"
    2926#include "psTest.h"
    3027#include "pslib.h"

  • trunk/psLib/test/dataIO/tst_psDB.c

    r3569 r3598  
    99 *  @author Aaron Culliney, MHPCC
    1010 *
    11  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-30 23:22:39 $
     11 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-31 23:01:46 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 */
    17 
    18 #define PS_ALLOW_MALLOC // for MySQL internals ...
    1917
    2018#include <unistd.h>
    … …  
    2624#include <string.h>
    2725
    28 #include "psDB.h"
    2926#include "psTest.h"
    3027#include "pslib.h"

  • trunk/psLib/test/dataIO/tst_psLookupTable_01.c

    r3330 r3598  
    1212*  @author  Ross Harman, MHPCC
    1313*
    14 *  @version $Revision: 1.12 $  $Name: not supported by cvs2svn $
    15 *  @date  $Date: 2005-02-25 19:55:01 $
     14*  @version $Revision: 1.13 $  $Name: not supported by cvs2svn $
     15*  @date  $Date: 2005-03-31 23:01:46 $
    1616*
    1717*  Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii
    … …  
    168168
    169169    // Allocate lookup table with valid parameters
    170     table1 = psLookupTableAlloc("verified/tableF32.dat", tableF32_validFrom, tableF32_validTo);
     170    table1 = psLookupTableAlloc("tableF32.dat", tableF32_validFrom, tableF32_validTo);
    171171    if(table1 == NULL) {
    172172        psError(PS_ERR_UNKNOWN,true,"Null lookup table generated from valid parameters");
    173173        return 1;
    174174    }
    175     if(strcmp(table1->fileName,"verified/tableF32.dat") != 0) {
     175    if(strcmp(table1->fileName,"tableF32.dat") != 0) {
    176176        psError(PS_ERR_UNKNOWN,true,"File name not properly stored in psLookupTable structure.");
    177177        return 2;
    … …  
    201201
    202202    // Allocate table using table with psU8 index and valid types and index-values
    203     table1 = psLookupTableAlloc("verified/tableU8.dat", 0, 100.5);
     203    table1 = psLookupTableAlloc("tableU8.dat", 0, 100.5);
    204204    if(table1 == NULL) {
    205205        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableU8.dat");
    … …  
    274274
    275275    // Allocate table using table with psS32 index and valid types and index-values
    276     table1 = psLookupTableAlloc("verified/tableS32.dat", -110, 1000.5);
     276    table1 = psLookupTableAlloc("tableS32.dat", -110, 1000.5);
    277277    if(table1 == NULL) {
    278278        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableS32.dat");
    … …  
    299299
    300300    // Allocate table using table with psF32 index and valid types and index-values
    301     table1 = psLookupTableAlloc("verified/tableF32.dat", -1100, 5500.5);
     301    table1 = psLookupTableAlloc("tableF32.dat", -1100, 5500.5);
    302302    if(table1 == NULL) {
    303303        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableF32.dat");
    … …  
    325325    // Allocate table using psU8 index but with unsorted rows and verify the list is
    326326    // sorted properly after being read
    327     table1 = psLookupTableAlloc("verified/table10.dat",1,2);
     327    table1 = psLookupTableAlloc("table10.dat",1,2);
    328328    if(table1 == NULL) {
    329329        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table10.dat");
    … …  
    397397
    398398    // Allocate table using table with invalid type in type row but valid index-values
    399     table1 = psLookupTableAlloc("verified/table2.dat",0,99.99);
     399    table1 = psLookupTableAlloc("table2.dat",0,99.99);
    400400    if(table1 == NULL) {
    401401        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table2.dat");
    … …  
    411411
    412412    // Allocate table using table with invalid value
    413     table1 = psLookupTableAlloc("verified/table3.dat",0,75.0);
     413    table1 = psLookupTableAlloc("table3.dat",0,75.0);
    414414    if(table1 == NULL) {
    415415        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table3.dat");
    … …  
    464464
    465465    // Interpolate values within the list and verify return values
    466     table1 = psLookupTableAlloc("verified/table10.dat",0,20);
     466    table1 = psLookupTableAlloc("table10.dat",0,20);
    467467    table1 = psLookupTableRead(table1);
    468468    for(psS32 i = 0; i < table1->numRows-1; i++ ) {
    … …  
    547547
    548548    // Interpolate values within the list and verify return values
    549     table1 = psLookupTableAlloc("verified/table10.dat",0,20);
     549    table1 = psLookupTableAlloc("table10.dat",0,20);
    550550    table1 = psLookupTableRead(table1);
    551551    statusVector = psVectorAlloc(table1->numCols,PS_TYPE_U32);

  • trunk/psLib/test/fileUtils/tst_psDB.c

    r3569 r3598  
    99 *  @author Aaron Culliney, MHPCC
    1010 *
    11  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-30 23:22:39 $
     11 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-31 23:01:46 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 */
    17 
    18 #define PS_ALLOW_MALLOC // for MySQL internals ...
    1917
    2018#include <unistd.h>
    … …  
    2624#include <string.h>
    2725
    28 #include "psDB.h"
    2926#include "psTest.h"
    3027#include "pslib.h"

  • trunk/psLib/test/fileUtils/tst_psLookupTable_01.c

    r3330 r3598  
    1212*  @author  Ross Harman, MHPCC
    1313*
    14 *  @version $Revision: 1.12 $  $Name: not supported by cvs2svn $
    15 *  @date  $Date: 2005-02-25 19:55:01 $
     14*  @version $Revision: 1.13 $  $Name: not supported by cvs2svn $
     15*  @date  $Date: 2005-03-31 23:01:46 $
    1616*
    1717*  Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii
    … …  
    168168
    169169    // Allocate lookup table with valid parameters
    170     table1 = psLookupTableAlloc("verified/tableF32.dat", tableF32_validFrom, tableF32_validTo);
     170    table1 = psLookupTableAlloc("tableF32.dat", tableF32_validFrom, tableF32_validTo);
    171171    if(table1 == NULL) {
    172172        psError(PS_ERR_UNKNOWN,true,"Null lookup table generated from valid parameters");
    173173        return 1;
    174174    }
    175     if(strcmp(table1->fileName,"verified/tableF32.dat") != 0) {
     175    if(strcmp(table1->fileName,"tableF32.dat") != 0) {
    176176        psError(PS_ERR_UNKNOWN,true,"File name not properly stored in psLookupTable structure.");
    177177        return 2;
    … …  
    201201
    202202    // Allocate table using table with psU8 index and valid types and index-values
    203     table1 = psLookupTableAlloc("verified/tableU8.dat", 0, 100.5);
     203    table1 = psLookupTableAlloc("tableU8.dat", 0, 100.5);
    204204    if(table1 == NULL) {
    205205        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableU8.dat");
    … …  
    274274
    275275    // Allocate table using table with psS32 index and valid types and index-values
    276     table1 = psLookupTableAlloc("verified/tableS32.dat", -110, 1000.5);
     276    table1 = psLookupTableAlloc("tableS32.dat", -110, 1000.5);
    277277    if(table1 == NULL) {
    278278        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableS32.dat");
    … …  
    299299
    300300    // Allocate table using table with psF32 index and valid types and index-values
    301     table1 = psLookupTableAlloc("verified/tableF32.dat", -1100, 5500.5);
     301    table1 = psLookupTableAlloc("tableF32.dat", -1100, 5500.5);
    302302    if(table1 == NULL) {
    303303        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with tableF32.dat");
    … …  
    325325    // Allocate table using psU8 index but with unsorted rows and verify the list is
    326326    // sorted properly after being read
    327     table1 = psLookupTableAlloc("verified/table10.dat",1,2);
     327    table1 = psLookupTableAlloc("table10.dat",1,2);
    328328    if(table1 == NULL) {
    329329        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table10.dat");
    … …  
    397397
    398398    // Allocate table using table with invalid type in type row but valid index-values
    399     table1 = psLookupTableAlloc("verified/table2.dat",0,99.99);
     399    table1 = psLookupTableAlloc("table2.dat",0,99.99);
    400400    if(table1 == NULL) {
    401401        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table2.dat");
    … …  
    411411
    412412    // Allocate table using table with invalid value
    413     table1 = psLookupTableAlloc("verified/table3.dat",0,75.0);
     413    table1 = psLookupTableAlloc("table3.dat",0,75.0);
    414414    if(table1 == NULL) {
    415415        psError(PS_ERR_UNKNOWN,true,"Unable to allocate table with table3.dat");
    … …  
    464464
    465465    // Interpolate values within the list and verify return values
    466     table1 = psLookupTableAlloc("verified/table10.dat",0,20);
     466    table1 = psLookupTableAlloc("table10.dat",0,20);
    467467    table1 = psLookupTableRead(table1);
    468468    for(psS32 i = 0; i < table1->numRows-1; i++ ) {
    … …  
    547547
    548548    // Interpolate values within the list and verify return values
    549     table1 = psLookupTableAlloc("verified/table10.dat",0,20);
     549    table1 = psLookupTableAlloc("table10.dat",0,20);
    550550    table1 = psLookupTableRead(table1);
    551551    statusVector = psVectorAlloc(table1->numCols,PS_TYPE_U32);

  • trunk/psLib/test/sysUtils/verified/tst_psConfigure.stderr

    r3130 r3598  
    66
    77<DATE><TIME>|<HOST>|I|psLibVersion00
    8     Current psLib version is: pslib-v1.4
     8    Current psLib version is: pslib-v1.5
    99
    1010---> TESTPOINT PASSED (psConfigure{Return current psLib version} | tst_psConfigure.c)
Note: See TracChangeset for help on using the changeset viewer.