Changeset 3497
- Timestamp:
- Mar 24, 2005, 12:36:16 PM (21 years ago)
- Location:
- trunk/psLib/src/astronomy
- Files:
-
- 2 edited
-
psAstrometry.c (modified) (4 diffs)
-
psAstrometry.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/astronomy/psAstrometry.c
r3264 r3497 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1. 59$ $Name: not supported by cvs2svn $11 * @date $Date: 2005-0 2-17 19:26:23$10 * @version $Revision: 1.60 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2005-03-24 22:36:16 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 26 26 #include "psConstants.h" 27 27 #include "psAstronomyErrors.h" 28 #include "psMatrix.h" 28 29 #include "psTrace.h" 29 30 #include "psLogMsg.h" … … 106 107 (X2 * -((E/B) / (F - ((C*E)/B)))) + 107 108 (Y2 * (1.0 / (F - ((C*E)/B)))); 109 110 XXX: Since thre is now a general psPlaneTransformInvertTmp() function, we 111 should rename this. 108 112 *****************************************************************************/ 109 113 static psPlaneTransform *invertPlaneTransform(psPlaneTransform *transform) … … 908 912 } 909 913 914 915 916 917 918 /***************************************************************************** 919 psPlaneTransformInvertHEY(out, in, region, nSamples): this is an earlier 920 version of this function which doesnot separate the code for building and 921 solving the matrixequations. 922 923 XXX: Delete this codeonce you know the other version works. 924 *****************************************************************************/ 925 psPlaneTransform *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 /***************************************************************************** 1064 multiplyCoeffs(trans1, trans2): Takes two 2-D polynomials as input and 1065 multiplies them. Basically, for each non-zero coeff in the trans1 coeff[][] 1066 array, you must multiply by all non-zero coeffs in trans2. 1067 1068 XXX: Inefficient in that the out polynomial is allocated every time. 1069 *****************************************************************************/ 1070 psDPolynomial2D *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 /***************************************************************************** 1102 psPlaneTransformCombineTmp(out, trans1, trans2) 1103 1104 XXX: Much room for optimization. Currently, we call the polyMultiply 1105 routine far too many times. 1106 *****************************************************************************/ 1107 psPlaneTransform *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 /***************************************************************************** 1242 psPlaneTranformFitTmp(trans, source, dest, nRejIter, sigmaClip) 1243 1244 XXX: What about nRejIter? Iterations? 1245 XXX: Use static vectors for internal data. 1246 *****************************************************************************/ 1247 bool 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 /***************************************************************************** 1362 psPlaneTransformInvertTmp(out, in, region, nSamples) 1363 1364 // XXX: Use static data structures. 1365 *****************************************************************************/ 1366 psPlaneTransform *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 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.3 4$ $Name: not supported by cvs2svn $11 * @date $Date: 2005-0 2-17 19:26:23$10 * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2005-03-24 22:36:16 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 534 534 ); 535 535 536 537 // XXX: These functions don't belong here. Will migrate to psCoords.c later. 538 // XXX: Doxygenate. 539 psPlaneTransform *psPlaneTransformInvertTmp( 540 psPlaneTransform *out, 541 const psPlaneTransform *in, 542 psRegion *region, 543 int nSamples 544 ); 545 546 // XXX: Doxygenate. 547 psPlaneTransform *psPlaneTransformCombineTmp( 548 psPlaneTransform *out, 549 const psPlaneTransform *trans1, 550 const psPlaneTransform *trans2 551 ); 552 553 // XXX: Doxygenate. 554 bool psPlaneTranformFitTmp( 555 psPlaneTransform *trans, 556 const psArray *source, 557 const psArray *dest, 558 int nRejIter, 559 float sigmaClip 560 ); 561 562 563 564 536 565 #endif
Note:
See TracChangeset
for help on using the changeset viewer.
