Changeset 6384
- Timestamp:
- Feb 8, 2006, 11:52:05 AM (20 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/astro/psCoord.c (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/astro/psCoord.c
r6383 r6384 10 10 * @author GLG, MHPCC 11 11 * 12 * @version $Revision: 1.10 8$ $Name: not supported by cvs2svn $13 * @date $Date: 2006-02-08 21: 34:53$12 * @version $Revision: 1.109 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2006-02-08 21:52:05 $ 14 14 * 15 15 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 58 58 "transform" is linear. 59 59 60 XXX: Since there is now a general psPlaneTransformInvert() function, we 61 should rename this. 60 XXX: This code no longer makes sense. The merge must be reviewed. 62 61 63 XXX: This code no longer makes sense. The merge must be reviewed. 62 XXX: below is the code using the standard matrix representation. note that 63 this inversion requires x->nX == 1, y->nY == 1 and x->nY <= 1, y->nX <= 1 64 64 *****************************************************************************/ 65 // XXX EAM : below is the code using the standard matrix representation.66 // note that this inversion requires x->nX == 1, y->nY == 1 and67 // x->nY <= 1, y->nX <= 168 65 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform) 69 66 { … … 177 174 } 178 175 179 // XXX: Must rewrite code and tests to use these functions.180 176 psPlane* psPlaneAlloc(void) 181 177 { … … 326 322 } 327 323 328 329 /******************************************************************************330 XXX: Private Function.331 332 XXX: Optimize this. We don't need a while loop here.333 334 piNormalize(): take an input angle in radians and convert it to the range 0:2*PI.335 *****************************************************************************/336 psF32 piNormalize(psF32 angle)337 {338 while (angle < FLT_EPSILON) {339 angle+=M_PI*2;340 }341 342 while (angle >= (M_PI*2)) {343 angle-=M_PI*2;344 }345 return(angle);346 }347 324 348 325 void projectionFree(psProjection *p) … … 455 432 } 456 433 457 psPlane* psProject(const psSphere* coord, 458 const psProjection* projection) 434 psPlane* psProject( 435 const psSphere* coord, 436 const psProjection* projection) 459 437 { 460 438 return(p_psProject(NULL, coord, projection)); … … 648 626 // 649 627 // XXX: rename, or verify, or recode, after the poly norder/nterm change. 628 650 629 psPlaneTransform *myPT = NULL; 651 630 if (out == NULL) { … … 791 770 792 771 XXX: What about nRejIter? Iterations? 793 XXX: Use static vectors for internal data.794 772 XXX: This code has problems with data that corresponds to a non-linear fit. 795 773 *****************************************************************************/ … … 805 783 PS_ASSERT_PTR_NON_NULL(dest, NULL); 806 784 785 psBool rc = true; 807 786 psS32 numCoords = PS_MIN(source->n, dest->n); 808 787 psS32 order = PS_MAX(trans->x->nX, trans->x->nY); … … 874 853 // Solution via LU Decomposition 875 854 // 876 // XXX: Check return codes877 //855 psVector *xSolution; 856 psVector *ySolution; 878 857 psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition 879 858 psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix 880 psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x 881 psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y 882 883 // 884 // Stuff coefficients into transformation 885 // 886 for (psS32 i = 0, ijIndex = 0; i < order; i++) { 887 for (psS32 j = 0; j < order - i; j++, ijIndex++) { 888 trans->x->coeff[i][j] = xSolution->data.F64[ijIndex]; 889 trans->y->coeff[i][j] = ySolution->data.F64[ijIndex]; 859 if (luMatrix == NULL) { 860 rc = false; 861 psError(PS_ERR_UNKNOWN, true, "psMatrixLUD() returned NULL. Returning FALSE.\n"); 862 } else { 863 xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x 864 ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y 865 866 if ((xSolution == NULL) || (ySolution == NULL)) { 867 rc = false; 868 psError(PS_ERR_UNKNOWN, true, "psMatrixLUSolve() returned NULL. Returning FALSE.\n"); 869 } else { 870 // 871 // Stuff coefficients into transformation 872 // 873 for (psS32 i = 0, ijIndex = 0; i < order; i++) { 874 for (psS32 j = 0; j < order - i; j++, ijIndex++) { 875 trans->x->coeff[i][j] = xSolution->data.F64[ijIndex]; 876 trans->y->coeff[i][j] = ySolution->data.F64[ijIndex]; 877 } 878 } 890 879 } 891 880 } … … 900 889 psFree(yVector); 901 890 902 return( true);891 return(rc); 903 892 } 904 893 … … 907 896 psPlaneTransformInvert(out, in, region, nSamples) 908 897 909 // XXX: Use static data structures.910 898 *****************************************************************************/ 911 psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out, 912 const psPlaneTransform *in, 913 psRegion region, 914 int nSamples) 899 psPlaneTransform *psPlaneTransformInvert( 900 psPlaneTransform *out, 901 const psPlaneTransform *in, 902 psRegion region, 903 int nSamples) 915 904 { 916 905 PS_ASSERT_PTR_NON_NULL(in, NULL); … … 1015 1004 if (rc == true) { 1016 1005 return(myPT); 1017 } 1018 1019 // XXX: Generate an error message, or warning message.1020 return(NULL);1006 } else { 1007 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not fit a transform to coordinates."); 1008 return(NULL); 1009 } 1021 1010 } 1022 1011
Note:
See TracChangeset
for help on using the changeset viewer.
