IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3497


Ignore:
Timestamp:
Mar 24, 2005, 12:36:16 PM (21 years ago)
Author:
gusciora
Message:

Added plane transform fit, invert, and combine functions.

These will eventually go in psCoord.c

Location:
trunk/psLib/src/astronomy
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/astronomy/psAstrometry.c

    r3264 r3497  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-02-17 19:26:23 $
     10 *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-03-24 22:36:16 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    2626#include "psConstants.h"
    2727#include "psAstronomyErrors.h"
     28#include "psMatrix.h"
    2829#include "psTrace.h"
    2930#include "psLogMsg.h"
     
    106107             (X2 * -((E/B) / (F - ((C*E)/B)))) +
    107108             (Y2 * (1.0 / (F - ((C*E)/B))));
     109 
     110XXX: Since thre is now a general psPlaneTransformInvertTmp() function, we
     111should rename this.
    108112 *****************************************************************************/
    109113static psPlaneTransform *invertPlaneTransform(psPlaneTransform *transform)
     
    908912}
    909913
     914
     915
     916
     917
     918/*****************************************************************************
     919psPlaneTransformInvertHEY(out, in, region, nSamples): this is an earlier
     920version of this function which doesnot separate the code for building and
     921solving the matrixequations.
     922 
     923XXX: Delete this codeonce you know the other version works.
     924 *****************************************************************************/
     925psPlaneTransform *psPlaneTransformInvertHEY(psPlaneTransform *out,
     926        const psPlaneTransform *in,
     927        psRegion *region,
     928        int nSamples)
     929{
     930    PS_PTR_CHECK_NULL(in, NULL);
     931    if (isProjectionLinear((psPlaneTransform *) in)) {
     932        return(invertPlaneTransform((psPlaneTransform *) in));
     933    }
     934    psPlaneTransform *myPT = NULL;
     935    psPlane *inCoord = psPlaneAlloc();
     936    psPlane *outCoord = psPlaneAlloc();
     937
     938    PS_PTR_CHECK_NULL(region, NULL);
     939    PS_INT_COMPARE(0, nSamples, NULL);
     940
     941    // XXX: Is this correct?
     942    psS32 order = PS_MAX(in->x->nX, in->x->nY);
     943
     944    //
     945    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
     946    //
     947    if (out == NULL) {
     948        myPT = psPlaneTransformAlloc(order, order);
     949    } else {
     950        if (!((out->x->nX == order) &&
     951                (out->x->nY == order) &&
     952                (out->y->nX == order) &&
     953                (out->y->nY == order))) {
     954            psFree(out);
     955            myPT = psPlaneTransformAlloc(order, order);
     956        } else {
     957            myPT = out;
     958        }
     959    }
     960    // XXX: Initialize myPT?
     961
     962    //
     963    // Create fake polynomial to use in evaluation
     964    //
     965    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
     966    for (int i = 0; i < order; i++) {
     967        for (int j = 0; j < order; j++) {
     968            fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
     969            fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
     970        }
     971    }
     972
     973    //
     974    // Create a grid of xin,yin --> xout,yout
     975    //
     976    psVector *xIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
     977    psVector *yIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
     978    psVector *xOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
     979    psVector *yOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
     980
     981    //
     982    // Initialize the grid of points
     983    //
     984    for (int yint = 0; yint < nSamples; yint++) {
     985        inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
     986        for (int xint = 0; xint < nSamples; xint++) {
     987            inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
     988
     989            (void)psPlaneTransformApply(outCoord, in, inCoord);
     990            xOut->data.F32[yint*nSamples + xint] = inCoord->x;
     991            yOut->data.F32[yint*nSamples + xint] = inCoord->y;
     992            xIn->data.F32[yint*nSamples + xint] = outCoord->x;
     993            yIn->data.F32[yint*nSamples + xint] = outCoord->y;
     994        }
     995    }
     996
     997    //
     998    // Initialise the matrix and vectors
     999    //
     1000    psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
     1001    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
     1002    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
     1003    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
     1004    for (psS32 i = 0; i < nCoeff; i++) {
     1005        for (psS32 j = 0; j < nCoeff; j++) {
     1006            matrix->data.F64[i][j] = 0.0;
     1007        }
     1008        xVector->data.F64[i] = 0.0;
     1009        yVector->data.F64[i] = 0.0;
     1010    }
     1011
     1012    //
     1013    // Iterate over the grid points
     1014    //
     1015    for (psS32 g = 0; g < nSamples*nSamples; g++) {
     1016        // Iterate over the polynomial coefficients, accumulating the matrix and vectors
     1017
     1018        for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1019            for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1020
     1021                fakePoly->mask[i][j] = 0;
     1022                psF64 ijPoly = psDPolynomial2DEval(fakePoly, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
     1023                fakePoly->mask[i][j] = 1;
     1024
     1025                for (psS32 m = 0, mnIndex = 0; m < order; m++) {
     1026                    for (psS32 n = 0; n < order - m; n++, mnIndex++) {
     1027                        fakePoly->mask[m][n] = 0;
     1028                        psF64 mnPoly = psDPolynomial2DEval(fakePoly, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
     1029                        fakePoly->mask[m][n] = 1;
     1030
     1031                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
     1032                    }
     1033                }
     1034
     1035                xVector->data.F64[ijIndex] += ijPoly * (psF64)xOut->data.F32[g];
     1036                yVector->data.F64[ijIndex] += ijPoly * (psF64)yOut->data.F32[g];
     1037            }
     1038        }
     1039    }
     1040
     1041    //
     1042    // Solution via LU Decomposition
     1043    //
     1044    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     1045    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
     1046    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
     1047    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
     1048
     1049    //
     1050    // Stuff coefficients into transformation
     1051    //
     1052    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1053        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1054            myPT->x->coeff[i][j] = xSolution->data.F64[ijIndex];
     1055            myPT->y->coeff[i][j] = ySolution->data.F64[ijIndex];
     1056        }
     1057    }
     1058
     1059    return(myPT);
     1060}
     1061
     1062
     1063/*****************************************************************************
     1064multiplyCoeffs(trans1, trans2): Takes two 2-D polynomials as input and
     1065multiplies them.  Basically, for each non-zero coeff in the trans1 coeff[][]
     1066array, you must multiply by all non-zero coeffs in trans2.
     1067 
     1068XXX: Inefficient in that the out polynomial is allocated every time.
     1069 *****************************************************************************/
     1070psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1,
     1071                                 psDPolynomial2D *trans2)
     1072{
     1073    psS32 orderX = (trans1->nX + trans2->nX) - 1;
     1074    psS32 orderY = (trans1->nX + trans2->nX) - 1;
     1075
     1076    psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
     1077    for (psS32 i = 0 ; i < out->nX; i++) {
     1078        for (psS32 j = 0 ; j < out->nY; j++) {
     1079            out->coeff[i][j] = 0.0;
     1080            out->mask[i][j] = 0;
     1081        }
     1082    }
     1083
     1084    for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) {
     1085        for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) {
     1086            if (0.0 != trans1->coeff[t1x][t1y]) {
     1087                for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
     1088                    for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
     1089                        out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
     1090                    }
     1091                }
     1092            }
     1093        }
     1094    }
     1095    return(out);
     1096}
     1097
     1098
     1099
     1100
     1101/*****************************************************************************
     1102psPlaneTransformCombineTmp(out, trans1, trans2)
     1103 
     1104XXX: Much room for optimization.  Currently, we call the polyMultiply
     1105routine far too many times.
     1106 *****************************************************************************/
     1107psPlaneTransform *psPlaneTransformCombineTmp(psPlaneTransform *out,
     1108        const psPlaneTransform *trans1,
     1109        const psPlaneTransform *trans2)
     1110{
     1111    PS_PTR_CHECK_NULL(trans1, NULL);
     1112    PS_PTR_CHECK_NULL(trans2, NULL);
     1113
     1114    //
     1115    // Determine the size of the new psPlaneTransform.
     1116    //
     1117    // PS_MAX(  Number of x terms in T2->x * number of x terms in T1->x,
     1118    //          Number of y terms in T2->x * number of x terms in T1->y,
     1119    psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX),
     1120                            (trans2->x->nY * trans1->y->nX));
     1121    psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY),
     1122                            (trans2->x->nY * trans1->y->nY));
     1123
     1124    psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX),
     1125                            (trans2->y->nY * trans1->y->nX));
     1126    psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY),
     1127                            (trans2->y->nY * trans1->y->nY));
     1128    psS32 orderX = PS_MAX(orderXnX, orderYnX);
     1129    psS32 orderY = PS_MAX(orderXnY, orderYnY);
     1130
     1131    //
     1132    // Allocate the new psPlaneTransform, if necessary.
     1133    //
     1134    psPlaneTransform *myPT = NULL;
     1135    if (out == NULL) {
     1136        myPT = psPlaneTransformAlloc(orderX, orderY);
     1137    } else {
     1138        if ((out->x->nX == orderX) && (out->x->nY == orderY) &&
     1139                (out->y->nX == orderX) && (out->y->nY == orderY)) {
     1140            myPT = out;
     1141        } else {
     1142            psFree(out);
     1143            myPT = psPlaneTransformAlloc(orderX, orderY);
     1144        }
     1145    }
     1146
     1147    //
     1148    // Initialize the new psPlaneTransform, if necessary.
     1149    //
     1150    for (psS32 i = 0 ; i < orderX ; i++) {
     1151        for (psS32 j = 0 ; j < orderY ; j++) {
     1152            myPT->x->coeff[i][j] = 0.0;
     1153            myPT->x->mask[i][j] = 0;
     1154            myPT->y->coeff[i][j] = 0.0;
     1155            myPT->y->mask[i][j] = 0;
     1156        }
     1157    }
     1158
     1159    //
     1160    // For each term (a * x^i * y^j) in trans2, we substitute the appropriate
     1161    // equation from trans1, and raise it to the appropriate power.  This is
     1162    // done via the multiplyDPoly2D().  The result is a polynomial (currPoly)
     1163    // and its coefficients are added into the myPT coeff matrix.
     1164    //
     1165    // XXX: This is horribly inefficient in that the trans1 polys are repeatedly
     1166    // multiplied against themselves.  This can easily be improved.
     1167    //
     1168    for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) {
     1169        for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) {
     1170            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
     1171            currPoly->coeff[0][0] = 1.0;
     1172            currPoly->mask[0][0] = 0;
     1173            psDPolynomial2D *newPoly = NULL;
     1174
     1175            if (trans2->x->mask[t2x][t2y] == 0) {
     1176
     1177                // Must raise trans1->y to the t2y-power.
     1178                for (psS32 c = 0 ; c < t2y; c++) {
     1179                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
     1180                    psFree(currPoly);
     1181                    currPoly = newPoly;
     1182                }
     1183
     1184                // Must raise trans1->x to the t2x-power.
     1185                for (psS32 c = 0 ; c < t2x; c++) {
     1186                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
     1187                    psFree(currPoly);
     1188                    currPoly = newPoly;
     1189                }
     1190
     1191                // Set the appropriate coeffs in myPT->x
     1192                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
     1193                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
     1194                        myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
     1195                    }
     1196                }
     1197            }
     1198            psFree(currPoly);
     1199        }
     1200    }
     1201
     1202
     1203
     1204    for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) {
     1205        for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) {
     1206            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
     1207            currPoly->coeff[0][0] = 1.0;
     1208            currPoly->mask[0][0] = 0;
     1209            psDPolynomial2D *newPoly = NULL;
     1210
     1211            if (trans2->y->mask[t2x][t2y] == 0) {
     1212
     1213                // Must raise trans1->y to the t2y-power.
     1214                for (psS32 c = 0 ; c < t2y; c++) {
     1215                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
     1216                    psFree(currPoly);
     1217                    currPoly = newPoly;
     1218                }
     1219
     1220                // Must raise trans1->x to the t2x-power.
     1221                for (psS32 c = 0 ; c < t2x; c++) {
     1222                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
     1223                    psFree(currPoly);
     1224                    currPoly = newPoly;
     1225                }
     1226
     1227                // Set the appropriate coeffs in myPT->x
     1228                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
     1229                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
     1230                        myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
     1231                    }
     1232                }
     1233            }
     1234            psFree(currPoly);
     1235        }
     1236    }
     1237
     1238    return(myPT);
     1239}
     1240
     1241/*****************************************************************************
     1242psPlaneTranformFitTmp(trans, source, dest, nRejIter, sigmaClip)
     1243 
     1244XXX: What about nRejIter?  Iterations?
     1245XXX: Use static vectors for internal data.
     1246 *****************************************************************************/
     1247bool psPlaneTranformFitTmp(psPlaneTransform *trans,
     1248                           const psArray *source,
     1249                           const psArray *dest,
     1250                           int nRejIter,
     1251                           float sigmaClip)
     1252{
     1253    PS_PTR_CHECK_NULL(trans, NULL);
     1254    PS_PTR_CHECK_NULL(source, NULL);
     1255    PS_PTR_CHECK_NULL(dest, NULL);
     1256
     1257    // Ensure that the input transformation is symmetrical.
     1258    if ((trans->x->nX != trans->x->nY) ||
     1259            (trans->y->nX != trans->y->nY) ||
     1260            (trans->x->nX != trans->y->nX)) {
     1261        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
     1262    }
     1263
     1264    psS32 numCoords = PS_MIN(source->n, dest->n);
     1265    // This is not really necessary because of above conditionals.
     1266    psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
     1267
     1268    //
     1269    // Create fake polynomial to use in evaluation
     1270    //
     1271    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
     1272    for (int i = 0; i < order; i++) {
     1273        for (int j = 0; j < order; j++) {
     1274            fakePoly->coeff[i][j] = 1.0;
     1275            fakePoly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
     1276        }
     1277    }
     1278
     1279    //
     1280    // Initialize the matrix and vectors
     1281    //
     1282    psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
     1283    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
     1284    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
     1285    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
     1286    for (psS32 i = 0; i < nCoeff; i++) {
     1287        for (psS32 j = 0; j < nCoeff; j++) {
     1288            matrix->data.F64[i][j] = 0.0;
     1289        }
     1290        xVector->data.F64[i] = 0.0;
     1291        yVector->data.F64[i] = 0.0;
     1292    }
     1293
     1294    //
     1295    // Iterate over the grid points
     1296    //
     1297    for (psS32 g = 0; g < numCoords; g++) {
     1298        // Iterate over the polynomial coefficients, accumulating the matrix and vectors
     1299
     1300        for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1301            for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1302                fakePoly->mask[i][j] = 0;
     1303                psF64 xIn = ((psPlane *) source->data[g])->x;
     1304                psF64 yIn = ((psPlane *) source->data[g])->y;
     1305                psF64 xOut = ((psPlane *) dest->data[g])->x;
     1306                psF64 yOut = ((psPlane *) dest->data[g])->y;
     1307                psF64 ijPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
     1308                fakePoly->mask[i][j] = 1;
     1309
     1310                for (psS32 m = 0, mnIndex = 0; m < order; m++) {
     1311                    for (psS32 n = 0; n < order - m; n++, mnIndex++) {
     1312                        fakePoly->mask[m][n] = 0;
     1313                        psF64 mnPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
     1314                        fakePoly->mask[m][n] = 1;
     1315
     1316                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
     1317                    }
     1318                }
     1319
     1320                xVector->data.F64[ijIndex] += ijPoly * xOut;
     1321                yVector->data.F64[ijIndex] += ijPoly * yOut;
     1322            }
     1323        }
     1324    }
     1325
     1326    //
     1327    // Solution via LU Decomposition
     1328    //
     1329    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     1330    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
     1331    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
     1332    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
     1333
     1334    //
     1335    // XXX: Should check the output of the matrix routines and return false if bad.
     1336    //
     1337
     1338    //
     1339    // Stuff coefficients into transformation
     1340    //
     1341    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
     1342        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
     1343            trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
     1344            trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
     1345        }
     1346    }
     1347
     1348    psFree(fakePoly);
     1349    psFree(permutation);
     1350    psFree(luMatrix);
     1351    psFree(xSolution);
     1352    psFree(ySolution);
     1353    psFree(matrix);
     1354    psFree(xVector);
     1355    psFree(yVector);
     1356
     1357    return(true);
     1358}
     1359
     1360
     1361/*****************************************************************************
     1362psPlaneTransformInvertTmp(out, in, region, nSamples)
     1363 
     1364// XXX: Use static data structures.
     1365 *****************************************************************************/
     1366psPlaneTransform *psPlaneTransformInvertTmp(psPlaneTransform *out,
     1367        const psPlaneTransform *in,
     1368        psRegion *region,
     1369        int nSamples)
     1370{
     1371    PS_PTR_CHECK_NULL(in, NULL);
     1372    //
     1373    // If the transform is linear, then invert it exactly and return.
     1374    //
     1375    if (isProjectionLinear((psPlaneTransform *) in)) {
     1376        return(invertPlaneTransform((psPlaneTransform *) in));
     1377    }
     1378    PS_PTR_CHECK_NULL(region, NULL);
     1379    PS_INT_COMPARE(0, nSamples, NULL);
     1380
     1381    // Ensure that the input transformation is symmetrical.
     1382    if ((in->x->nX != in->x->nY) ||
     1383            (in->y->nX != in->y->nY) ||
     1384            (in->x->nX != in->y->nX)) {
     1385        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
     1386    }
     1387    psS32 order = PS_MAX(in->x->nX, in->x->nY);
     1388
     1389    psPlaneTransform *myPT = NULL;
     1390    psPlane *inCoord = psPlaneAlloc();
     1391    psPlane *outCoord = psPlaneAlloc();
     1392
     1393    //
     1394    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
     1395    //
     1396    if (out == NULL) {
     1397        myPT = psPlaneTransformAlloc(order, order);
     1398    } else {
     1399        if ((out->x->nX == order) && (out->x->nY == order) &&
     1400                (out->y->nX == order) && (out->y->nY == order)) {
     1401            myPT = out;
     1402        } else {
     1403            psFree(out);
     1404            myPT = psPlaneTransformAlloc(order, order);
     1405        }
     1406    }
     1407
     1408    //
     1409    // Copy the input transform to myPT.
     1410    //
     1411    for (psS32 i = 0 ; i < in->x->nX ; i++) {
     1412        for (psS32 j = 0 ; j < in->x->nY ; j++) {
     1413            myPT->x->coeff[i][j] = in->x->coeff[i][j];
     1414        }
     1415    }
     1416    for (psS32 i = 0 ; i < in->y->nX ; i++) {
     1417        for (psS32 j = 0 ; j < in->y->nY ; j++) {
     1418            myPT->y->coeff[i][j] = in->y->coeff[i][j];
     1419        }
     1420    }
     1421
     1422    //
     1423    // Create a grid of xin,yin --> xout,yout
     1424    //
     1425    psArray *inData = psArrayAlloc(nSamples * nSamples);
     1426    psArray *outData = psArrayAlloc(nSamples * nSamples);
     1427    for (psS32 i = 0 ; i < inData->n; i++) {
     1428        inData->data[i] = (psPtr *) psPlaneAlloc();
     1429        outData->data[i] = (psPtr *) psPlaneAlloc();
     1430    }
     1431
     1432    //
     1433    // Initialize the grid.
     1434    //
     1435    psS32 cnt = 0;
     1436    for (int yint = 0; yint < nSamples; yint++) {
     1437        inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
     1438        for (int xint = 0; xint < nSamples; xint++) {
     1439            inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
     1440            (void)psPlaneTransformApply(outCoord, in, inCoord);
     1441
     1442            ((psPlane *) outData->data[cnt])->x = inCoord->x;
     1443            ((psPlane *) outData->data[cnt])->y = inCoord->y;
     1444            ((psPlane *) inData->data[cnt])->x = outCoord->x;
     1445            ((psPlane *) inData->data[cnt])->y = outCoord->y;
     1446
     1447            cnt++;
     1448        }
     1449    }
     1450    bool rc = psPlaneTranformFitTmp(myPT, inData, outData, 10, 100.0);
     1451
     1452    psFree(inCoord);
     1453    psFree(outCoord);
     1454    psFree(inData);
     1455    psFree(outData);
     1456
     1457    if (rc == true) {
     1458        return(myPT);
     1459    }
     1460
     1461    // XXX: Generate an error message, or warning message.
     1462    return(NULL);
     1463}
  • trunk/psLib/src/astronomy/psAstrometry.h

    r3264 r3497  
    88*  @author GLG, MHPCC
    99*
    10 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2005-02-17 19:26:23 $
     10*  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2005-03-24 22:36:16 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    534534);
    535535
     536
     537// XXX: These functions don't belong here.  Will migrate to psCoords.c later.
     538// XXX: Doxygenate.
     539psPlaneTransform *psPlaneTransformInvertTmp(
     540    psPlaneTransform *out,
     541    const psPlaneTransform *in,
     542    psRegion *region,
     543    int nSamples
     544);
     545
     546// XXX: Doxygenate.
     547psPlaneTransform *psPlaneTransformCombineTmp(
     548    psPlaneTransform *out,
     549    const psPlaneTransform *trans1,
     550    const psPlaneTransform *trans2
     551);
     552
     553// XXX: Doxygenate.
     554bool psPlaneTranformFitTmp(
     555    psPlaneTransform *trans,
     556    const psArray *source,
     557    const psArray *dest,
     558    int nRejIter,
     559    float sigmaClip
     560);
     561
     562
     563
     564
    536565#endif
Note: See TracChangeset for help on using the changeset viewer.