IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

cosmetic tweaks.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.